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ABSTRACT 


Tracking a maneuvering acoustic source using multipath (MP) arrivals in an 
inhomogeneous (IH) ocean medium is investigated. Errors introduced by a horizon- 
tally stratified sound speed profile are quantitatively evaluated. A new method of 
converting MP time difference of arrivals to depth and range which accounts for the 
IH effect is developed and evaluated. A 3-D target tracker previously used in MP 
tracking is modified in order to remove estimation biases and improve computa- 
tional efficiency. Tracking performance is demonstrated using extensive simulation 


and shown to be greatly improved. 
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I. INTRODUCTION 


A. BACKGROUND - PASSIVE TARGET TRACKING 
Passive tracking is based on reception and analysis of signals emitted from 
a target at an unknown time and location. Range measurement to the target, 
usually based on knowledge of the travel time and speed of the signal, is not 
possible. The lack of range information makes the tracking of targets, which are 
free to change their range, a challenge. Our interest in passive target tracking is 
primarily for acoustic tracking of targets in the ocean. However, the problem of 
passive localization is encountered in fields like radio astronomy and seismology as 
well as in passive sonar. Several passive tracking techniques have evolved over the 
years; some of these will be briefly reviewed here. 
1. T.M.A - Target Motion Analysis 
Target motion analysis is a range estimation technique devised for World 
War II submarines. Typical targets were surface ships and the measurements 
were primarily sonar bearing. Occasional course and speed estimates derived from 
periscope peeks were used as well, but mainly for initialization. The method is 
based on hypothesizing constant course and speed targets. The range of the tar- 
get whose trajectory woe have produced bearing measurements which best fit 
the actually observed bearing data. is selected as the estimate. The observability 
problem is partially solved by maneuvering the observing platform while assuming 
the target maintains its course and speed. Own ship maneuvering is a lengthy 
operation and the assumption is that the target maintains course and speed only 


of limited validity. The important points to note about TMA are: 


20 


e No direct or indirect measurement of range. 
e Range estimates limited to nonmaneuvering targets. 
e Lengthy own ship maneuver requirement. 
e Lengthy solution development time. 
e Very accurate bearing measurements needed. 
2. Wavefront Curvature 

Wavefront curvature is an indirect range measurement technique. It is 
- based on the assumption that the acoustic waves propagate in spherical wavefronts. 
Signals from an array of at least three spatially separated receivers _ delayed in 
time to affect focusing . The amount of delay, which is the travel time difference 
of arrival (TDOA) required for a given array length, is then translated into target 
range. The performance depends on the length of the array baseline, the precision 
of both the sensor location along the array and the time delay measurements. 

A number of implementations exist with array baselines varying from a 
few tens of meters in integrated systems, to hundreds of meters in large multi- 
platform distributed configurations. The range measurement, which is indepen- 
dent, can however be further combined with other measurements like target bear- 
ing or Doppler in order to reduce the ranging error. The important features of this 


method are: 
e Truly independent measurement of range. 
e Requirement for multiple sensors precisely located over a large baseline. 
e Applicability to maneuvering targets. 


3. Multipath Tracking 


Multipath (MP) tracking is similar in principle to wavefront curvature. 


Here travel time differences to a single receiver via different paths, are measured 
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using autocorrelation techniques. The paths most commonly used are the direct 
path and those reflected from the surface and bottom of the ocean. The reflected 
paths can be viewed as paths to imaginary virtual receivers in mirror image lo- 
cations, one above the water and the other below the bottom. In this sense, the 
single receiver acts like an array of three receivers (see Section C.2). 

As in the wavefront methods, the travel time differences depend on the 
target’s depth and range and on the array baseline, here determined by the water 
column. If propagation is assumed to be along straight lines, the inverse dependence 
of target position on travel time differences is relatively simple. It provides a 
mapping of time delays to range and depth. The range and depth observations 
thus obtained can again be combined with other available target measurements in 
order to reduce the error. The main features of multipath tracking are therefore: 

e Independent measurement of range and depth. 
e Compact single receiver configuration. 

e Applicability to maneuvering targets. 

e Dependence on ocean propagation conditions. 

The MP method has enjoyed a growing amount of attention in recent years 
due to the advantages stated above and the continued improvement in time delay 
estimation. Compensating for the dependence of the method on ocean conditions 


is the main subject of this work. 


B. THE BASIC MULTIPATH METHOD 
1. The Travel Time Differences 
The acoustic medium of the ocean is confined between the distinct surface 
and bottom boundaries. The medium thus provides multiple reflective paths along 


which sound can travel between two given points. The sound travel time along 


om 


those different paths will differ by an amount dependent in a known manner on 


the geometry of the scenario. 


Consider the direct and single reflection paths shown in Fig. 1.1, under 


the following assumptions: 





TO 


Fig. 1.1 Multipath propagation. 








ee ee 


A. The sound propagates along straight lines. 


B. The bottom and surface are ideal rigid and pressure release boundaries, re- 
spectively. 


C. The travel time differences can be analyzed, resolved, and associated with the 
corresponding paths. 


The following relations then exist between the target depth D and range R position, 


measured relative to the observer, and the travel time differences 7 and 72: 


n =T,-T) = 5: [(R? + D? +4D3 — 4D)‘ — (1-14) 


and 


T = Tz, —To = 2 (FR MD? 44(Ds — Dg)? 4D, — Do)D)” - p|(1.16) 


Where 7o,71,7> are the travel times of the direct, surface, and bottom paths, 
respectively; D,, is the water depth ( surface to bottom); D, is the observer depth 
(surface to observer); p is the slant range given by (R? + D?)1/? and C is the speed 
of sound. These relation can be inverted to express the slant range, the depth and 


the range as functions of the travel time differences. The results are 


ATO S= Dy) 2 DDS —D,)") — C?T? Dy + T?(Dw — Da) 
i SO Denno. To Dal 


- |C?Ty + 2CT2p — 4(Dwy — Do)*| (1.25) 


(1.2a) 


1 


om AD neDs | 


and 


1/2 


R = (p* — D?) (lec) 


Eq. (1.1) and Eq. (1.2), referred to in this work as the direct and inverse functions 
respectively, were derived by Hassab in Ref. 1. 
2. Time Delay Estimation 
The signal y’(t) at the receiver is modeled as a sum of replicas of the 
original signal z(t), each multiplied by the path gain a;, and time shifted by the 
path travel time T;, that is, 
N 


y(t) =} a;-2(t-T;) es) 


2=0 


ao 


The signal y(t), which is the signal y’(t) time shifted by the travel time of the first 
arrival To, that is, 
y(t) = y'(t — To) (1.4) 


can be expressed in terms of the TDOAs 7; as 


N N 
y(t) = > aiztt —T) —T;)= > aia(t — 7;) (1.5) 


where 7; = 7; — To. If the source of the target’s acoustic emission is propeller 
cavitation noise, as is the case for broadband MP, then the original signal z(t) is 
a broadband random process. The autocorrelation of its received version y(t) will 
exhibit peaks at time lags equal to 7;. The autocorrelation can be easily computed 
at discrete time lags using a sampled version of the received signal. As long as 
the ACF peaks are resolvable, they can be interpolated between the discrete lags 
providing the required continuous time delay measurement. The measured delays 
can be further smoothed using a time delay tracker. 

In the narrowband case, low frequency tonals emitted by target machinery 
can be received at long range. Since the autocorrelation function of such sinusoids is 
periodic, the ACF approach is impractical. A different approach has been applied 
to this situation which is based on standing waves generated by the tonals in 
the water column. The sound field intensity is sensed by a number of receivers 
placed along a vertical array. A technique called field matching is used to select a 
target position which will give rise to a sound field that best matches the measured 
intensity. The effort in this work concentrates on the broadband case. 

3. Target Tracking | 
The MP depth and range measurements are contaminated with noise. 


Sources of the noise include acoustical noise which distorts the autocorrelation 
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of the received signal, errors and fluctuations of the MP conditions in the ocean, 
and errors in the receiver and signal processor. To estimate the actual position 
from the noisy measurements, a target tracker is employed. The tracker fits the 
measurements with an estimate that, on the average, minimizes the squared error 
between the actual and estimated position. Such a tracker, also referred to as 
an estimator, observer, or filter, may also use other available measurements like 
bearing and Doppler. 

Many trackers of the type described above are available. These range 
from simple averaging trackers to large banks of Kalman filters. The trackers differ 
primarily in the measurements they use and the complexity of the assumed target 
model. 

A typical complete MP measurement and tracking system is shown in 
Fig. 1.2. Estimated TDOAs are converted to depth and range by the prefilter. 
Depth and range are then combined with bearing to form a 3-D target position 
measurement which is filtered by the target tracker. 


Current MP tracking methodology is lacking in the following areas: 


e The methods assume straight line propagation which is true only in a homoge- 
neous medium. This assumption is largely in error for the realistic inhomoge- 
neous (IH) case. This causes large errors in the transformation of time delays 
to depth and range. 


e The assumption that path and delays are easily associable is wrong for many 
practical cases. A significant amount of detailed oceanographic data is re- 
quired in order to identify the multiple paths. This information is not typically 
available in an actual target tracking situation. 


e Receiver time delay resolution is not perfect but limited by the bandwidth 
of the signal and the receiver. This affects the performance whenever the 
geometry yields two or more multipaths with similar travel times. Three such 
typical geometries are : 
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Fig. 1.2. MP tracking system. 


o Long range, where all paths tend to be of same length. 


o Target or observer close to the surface or bottom where the direct path 
and the reflected path have similar length. 


o Target and observer at symmetrically opposing depths such that the bot- 
tom and surface paths have similar travel times. 


e The lack of exact knowledge of bottom depth and structure limits both the 
time delay measurement and the exact translation to depth and range. 


e Other effects including slanted bottom, and the anisotropic nature of the bot- 
tom which may have different slopes at different directions. 


Our work seeks to address some of these effects, in particular the first 


three. 
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C. REVIEW OF PAST WORK 
1. MP Target Tracking 
a. Hassab 
In 1976 Hassab published an original work titled “Passive tracking of 
moving source by a single observer in shallow water” [Ref. 1] which was the starting 
point for the research in the area. In this work, the use of measured travel time 
differences as inputs to a target tracking filter was established. The tracker was 
realized in cartesian coordinates using an extended Kalman filter and assuming a 
nonmaneuvering target. Earlier work by Hassab [Ref.2] dealt with the very process 
of measuring the travel time differences using a technique based on the cepstrum. 
b. Singer 
In the early T.M.A algorithm the development of a range solution was 
halted when a target maneuver was detected (a “Zig Zag” in World War II jargon). 
The maneuvers were detected when constant tracking errors indicated a mismatch 
between the actual target and its model. With the application of Kalman filters 
to target tracking, the maneuvering command uncertainty was represented as a 
white gaussian process noise. Singer [3] improved the model by coloring (low pass 
filtering) the command noise to correspond to the expected maneuver dynamics. 
c. Moose 
In the 1970’s multiple model (multiple hypothesis) methods were in- 
troduced [4] to address the maneuvering problem. Here a number of target models 
based on different hypothesized maneuvers are computed in parallel, and combined 
in a Bayesian manner to form the overall estimate. 
This approach was applied to both airborne and underwater targets by 
Moose [5,11] and his colleagues McCabe [6], Gholson [7], Van Landingham [8], Dai- 


ley [9], and Caputi {10]. Three main issues related to MP tracking were addressed 
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in this work, which is ongoing; They are the selection of a particular coordinate 
system, the account for target maneuver, and the bias associated with the nonlin- 
ear transformation of TDOAs to depth and range. Performance evaluations were 
carried out in Moose’s work but both the simulation and tracking were done with- 
out taking account for the effects of the IH acoustic medium and the realistic delay 
estimation process. 

2. Time Delay Estimation and Source Localization. 

The area of time delay estimation has received enormous attention in the 
last twenty years. Among the contributors are Schultheiss [13], Ianniello [14], Wien- 
stien {15], Friedlander [16, 17], and others. Work in this area covered the topics 
of estimation instrumentation, theoretical and experimental error bounds, estima- 
tion resolution, continuous delay reconstruction from sampled time sequences, and 
more. A very good summary is presented in Refs. 15 and 16. 

Use of TDOA for source localization was investigated by Schultheiss for a 
variety of receiving array configurations. In his work [18] lower bounds on achiev- 
able error using bandlimited signals is developed using the Cramer Rao lower bound 
(CRLB) and Ziv Zakai lower bound (ZZLB). In the Ph.D. dissertation by his stu- 
dent Hamilton [19] the lower bound for the MP localization error variance is also 
developed. Hamilton also establishes the relation between the wavefront curvature 
and the multipath ranging using reflected mirror images of the receiver. 

3. Acoustic Propagation and Modeling 

The basic form of sound propagation in an inhomogeneous medium where 
the speed of sound is a function of depth has long been known. The development 
of the associated approximations and the resulting ray acoustics are presented, for 


example, in Ref. 20 by Ziomek. 


34 


Numerous computer models have been developed to compute sound prop- 
agation in the ocean based on ray acoustics, starting from analog computer sim- 
ulations which were used for ray tracing, and proceeding up to very large and 
sophisticated numerical models which account for many of the special effects of 
the ocean medium. A classic model in this category is the Generic Sonar Model 
developed by Weinberg [21]. 

Only a few models address the problem of finding the rays traveling be- 
tween two given end points, the eigenray problem. A recent model of this kind is 
SMART (SMall Acoustic Ray Tracer) developed by Novick [22]. 

Research on the nature of reflection from the bottom and the surface of 
the ocean also dates back to World War II. Contributors in this field include, for 
example, Clay and Medwin [23]. Recent work in this area reveals the interrela- 
tionship between the various oceanographic processes, for example, internal waves, 
surface interaction of the ocean and the atmosphere, and shear waves induced in 
the bottom of the ocean by a sound field. The limitation of simplified partial, 


lumped models for these phenomena is becoming apparent. 


D. PROBLEM STATEMENT 


This thesis research seeks to address the following issues: 


e To account for the complex effects of the inhomogeneous (vertically stratified ) 
ocean medium on broadband MP tracking. 


e To investigate the impact and account for the effects of a realistic receiver and 
time delay estimator on the overall tracking process. 


In addition our work includes: 


e Modification of the state of the art target tracking algorithm to decrease some 
of its estimation bias. 
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e Evaluation of the complete tracking algorithm using a realistic environment, 
and target simulation considering all main sources of estimation errors. 


The main contributions and key discussions of each of the foregoing issues are given 
below. 
1. The Inhomogeneous Ocean Medium 

The real ocean 1s acoustically inhomogeneous , in that among other effects, 
the speed of sound varies with depth. This variation significantly effects the MP 
propagation and travel time as shown in Fig. 1.3 where the receiver (Rr) is at 
range = 0. Note that the direct path between the source and the receiver is com- 
pletely eliminated due to the ray bending in Fig. 1.3b (an exhaustive search for the 
eigen rays was conducted for this plot and the resultant eigenrays are plotted). The 
effects, believed to be analyzed quantitatively here for the first time, are shown to 
render the assumption of straight line propagation to be of limited practical use. 
Accounting for this effect is difficult since the inverse function, which transforms 
travel time differences to depth and range is not readily computable by existing 
acoustic models. Devising an inversion method to account for the IH effects on the 
MP tracking 1s a main goal of this research. 

The acoustic energy attenuates as it propagates through the water, and 
many factors contribute to this loss including spreading, absorption, and reflec- 
tion. While these loses and their impact on the time delay estimation noise are 
relatively well understood, they were nevertheless never considered in past simula- 
tions. Inclusion of the reflection and spreading effects to increase the reliability of 
the simulation is another part of the effort to better account for the effects of the 


Medium. 
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Fig. 1.3. IH MP propagation, a. with direct path, b. without direct path. 
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2. Realistic Receiver and Delay Estimation 


a. Limited Resolution 
The finite bandwidth of the target signal and receiver limit the time 
delay resolution to a few tenths of a millisecond. A significant amount of this 
research was dedicated to prediction and improvement of this limitation but less 
work was aimed at investigating the details of its impact on MP tracking. The 
effect of this limited resolution on the MP tracker and the means to overcome it 
are another subject addressed by this research. 
b. Nonidentifiable Paths 
The polarity of the ACF can support association of the time delays 
with their corresponding paths when the MP structure is known to be simple. 
However, when the structure is complex and/or unknown, such as in cases where 
there is a lack of direct path due to ray bending, the simple association is not 
possible. The ability to associate delay with paths, assumed in previous work, is 
not assumed here. 
c. Nonnegative Noise 
The even symmetry of the autocorrelation function does not enable 
the measurement of time delay polarity. In the multipath situation, there is a first 
arrival followed by lagging replicas. This justifies the use of the positive time delays 
only. The implication is that the delay noise is not normally distributed as was 
assumed in the past. A distribution which is zero for negative delay is a better 
description and it leads to a nonzero delay estimation bias error. 
3. Tracker Modification 
In the maneuvering target tracker radar by Moose, the conditional 
mean of some hypotheses commands is used as the estimate. This estimate is 


biased if the hypotheses commands are not symmetrically centered around the 
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actual value. Removing this bias is difficult since modification of the hypotheses 
disrupts a recursion used in the estimation process. Modification of the tracker to 
remove the bias without disrupting the recursion or introducing large estimation 
transients is one of the goals of the research. 
4. Evaluation 
Many elements are involved in MP tracking, namely 
e The target and its dynamics. 


e The medium and its multipath structure. 


The receiver and the time delay analyzer. 
e The conversion from time delays to position measurement. 
e The actual acquisition and tracking of a maneuvering target. 

Evaluation of the system performance is done in two steps. First, each and 
every component is evaluated separately. This is done under simulated conditions, 
which are determined for each component based on the analysis of the overall 
system operation. Then and only then, the integrated system is evaluated, using 


the results of the individual component evaluation as ‘reference data’. 


E. SCOPE AND OUTLINE 

This work approaches the MP tracking as a system problem. The remaining 
chapters are as follows. 

Chapter Two deals with the 3-D target tracker. A state-of-the-art tracker 
is described. An idealized receiver delay analyzer and homogeneous straight line 
propagation are assumed. The algorithm provides a good starting point in terms 
of its treatment of a 3-D maneuvering target. New interpretation and analysis is 
provided for certain aspects of the estimation mechanism, leading to a modification 


and improvement . A systems approach, not used in previous work, is applied to 
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evaluate the performance of the modified tracker. Some of the parameters and 
algorithms used for the simulation are not discussed in detail here, since they 
result from the removal of the above idealized assumptions which is explained in 
Chapter Three. 

Chapter Three addresses the main source of tracking error, namely the JH 
medium together with the effects of a realistic receiver and delay analyzer. The 
nature of the problem is described and the proposed solution is discussed in detail. 
The performance of the method is then evaluated and analyzed in detail. 

Chapter Four presents a detailed evaluation and analysis of the performance 
of the new inversion prefilter. 

Chapter Five resumes the overall See performance analysis, this time with- 
out idealized assumptions and using the new inversion algorithm and the improved 
tracker. 

The work is summarized in Chapter Six where conclusions and recommended 


extensions are presented. 
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Il. THREE DIMENSIONAL MANEUVERING TARGET TRACKER 


A. INTRODUCTION 

The three dimensional (3-D) target tracker developed in this chapter is based 
on a state-of-the-art design by Moose, and McCabe. The modified Kalman filter 
which forms the central portion of the tracker was recently investigated by Saez 
[Ref. 24], in conjunction with this study. In his work, Saez includes a detailed 
review of the tracker, results of which are briefly repeated here for completeness 
and to establish notation. 

The emphasis in this section is on improvements of the tracker design and a 
new overall systems approach to the performance evaluation. The modifications in- 
clude an adaptive, instead of a fixed, command hypothesis bank and an advancing 
smoother. Both modifications are intended to reduce tracking biases. The investi- 
gation also includes the use of a second instead of a third order target model. This 
was found to reduce computational load without sacrificing estimation accuracy. 

A realistic new model for the time delay estimation noise is used in evaluating 
the performance. The model incorporates propagation effects as well the effects 
of some inaccuracies in the time delay estimation. This enables a more realistic 
evaluation of the overall performance for the homogeneous case, which will be 
extended in Chapters Three through Five to the IH case. 

A range and bearing coordinate decoupling approximation was introduced in 
the original tracker by McCabe [Ref. 6] to linearize the model. An interesting 
interpretation of this procedure provides an explanation for errors in tracking non- 


maneuvering targets that occur at short ranges. 
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B. TARGET DESCRIPTION 
1. Horizontal Plane 
A second order physical model is used to describe target motion in the X 
and Y directions of the horizontal plane. Contrclling the motion is a command 
which forms the third degree of freedom in each axis. If Th(t) is the command 
thrust in Newtons (N) and D, is the drag in [N per m/sec] the differential equation 


for each axis which follows from Newton’s second law is 
Th(t) — D.z(t) = mz(t). (2.1) 
This equation can be rewritten as 


z(t) = —a (U(t) — z(t) (2.2) 


where a = D,/m [sec~"] is the reciprocal of the system’s damping time constant, 
and U = Th(t)/D, [m/sec] is the speed command, i.e., the speed at which the plat- 
form will move when steady state is reached. If it is assumed that X(0) = X(0) =0 


and that the command is a constant U then the solution is : 


z(t) = U(1 —e7**) (2.3) 


a(t) =U (t+ 2 


A (e7o* — .)) (2.4) 


2. Depth 
Because of the reduced dynamics expected for underwater targets in the 
depth channel, the target’s motion in depth was modeled as a first order system. 
The command Uy {m] can therefore be represented in terms of the steady state 


depth, and the differential equation along the depth axis is 
D(t) = aa(Ug — D(t)). (2.5) 
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The solution of Eq.(2.5) is given by 
De) = Dever sp Ube ul = haa) (2.6) 


3. Command Noise Augmentation 


Eq. (2.2) describes target motion in response to a known control U. How- 
ever the control (command) is not known to the observer, especially if the target 1s 
maneuvering. Singer partially accounts for this by adding a colored noise compo- 
nent w’ to the command. This = w’ is modeled by a first order recursive filter 


driven by a normally distributed white noise input w. 
Wi) = —Aww'(t) + w(t) (2.7) 


The lowpass model was chosen since it represents a maneuver which typically takes 
at least a few seconds to complete. The command is thus modeled in the ob- 
server by the sum of an assumed command U and the random noise w’, that is, 
Uosserver = U + w'(t). This description suggests estimation of the command U by 
the multi-model (MM) estimator described in Section C. 

In order to maintain a complete system state description, the coloring 
LP filter and its state w’ are combined into the target state equations. Both the 
command U and the the noise w are represented in units of the resulting steady 
state speed. The differential equations are then represented in discrete form by the 


following set of (3-D) Cartesian motion difference equations. 


Xn = OXn-1 + BL. + Vw, n-1 (2.8a) 
es = YnPY¥ n-1 + YU, + Vw, n—1 (2.85) 
D, = aa Dn-; + (1 _ aqg)U4q (Z23c) 
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where 7 
K=( 2 


Y=(y 9 wy)” 
and the components of the matrices ¢,I,W are detailed in Appendix A, and 
aq = e~%4?. Note that since the equations are linear and decoupled, a change 
in target course at a constant speed is modeled here as an acceleration along one 

axis and a deceleration along another. 

4. Cylindrical Coordinate Observations 
Most surveillance and tracking sensors produce measurements relative to 
their Vw.. position and orientation. The same 1s true in the MP case, where range 
and depth are indirectly deduced from the travel time differences across the vertical 
plane. Bearing is measured by beamforming; here sonar azimuthal beams are 
formed using travel time differences across the horizontal plane. Thus, the natural 
measurements are in the cylindrical coordinate system of range, depth and bearing. 
These are relative to the observing sensor which is therefore set at the origin. The 
detailed set of variables used to describe the positions of the observer and the target 


in the ocean is shown in Fig. 2.1 and defined in Table 2.1. 
TABLE 2.1 


CYLINDRICAL COORDINATES VARIABLES 


Do Observer depth measured from surface. 

Dw Depth of the water. 

D Depth of target relative to observer. 

R Horizontal range of target relative to observer. 
B Bearing (azimuthal target angle relative to north). 
B Bearing rate. 

U, Speed command in range direction. 

Ub Speed command in cross range direction. 

Ud Depth position command. 

U, Depth speed ag(Ug — D) 

p Slant range. 
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Fig. 2.1 Cylindrical coordinates. 


5. Linearization and Decoupling 


In order to use the cylindrical measurements as observations for the Carte- 
sian target model of Eq. (2.8), the measurements have to be transformed. This 
leads to a coupled, nonlinear observer. Two possible approaches have been applied 
to this problem in the past. One is based on conversion and linearization of the 
measurements, and leads to the extended Kalman filter, (EKF). This approach was 
applied to MP tracking by Hassab [1]. The other was introduced by Moose [5] and 
is based on a decoupling and linearization of the state equations, using discretized 
approximations. 

Investigation of the two approaches [7] indicated a preference for the second 
since taking account of target maneuvers in the EKF seemed to require a very large 
computational load. The second approach, developed in detail by McCabe in [6], 
is based on the following main steps. 

The Cartesian target state description is transformed to cylindrical coor- 


dinates using the transformation: 


R= (zr? +y2)'? (2.9a) 
B =waa-* G/z) (2.95) 
DvD. (2.9c) 


In the bearing channel, the range is assumed constant for the duration of 
the sampling period. The speed maneuver commands U, and U, along the X and 
Y directions are replaced by commands along the range and cross range (bearing) 
coordinates U, and U, respectively. 

Range and bearing coordinate coupling is thus reduced to a parametric 


relation. This requires the bearing channel matrices to be recomputed only when 
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the range changes by more then some given ratio and no more then once per 


iteration. The resulting state equations are 


Dn = agDn-1 + (1 — ag) Ud,n-1 (2.10a) 
R, = $-Rn-1 +T U0, n-1 + Vrw, (2.105) 
Bh = %Bn-1 +PU, n-1 + Vows (2.10c) 
where 
R=(RRw')* (2.11) 
and 
B=(BBw,). (2e12) 


Each channel of the three decoupled channels (depth, range, bearing) is observed 
by a scalar z defined as 


Zdn = dn + VoIn (2.13a) 
ee Ue (2.136) 
Zpn = Bat Vin (2.13c) 


where vq, v,, and v, are the observation noises. The components ¢, I, WY, of 
the cross range matrices are range dependent and presented in Appendix A. The 
subtle terminology — preference of the term “cross range” to “bearing” — used by 
McCabe, is significant. It is indicative of the fact that while the bearing variable 
is used, a constant cross range rate motion is actually being modeled as the basic 
nonmaneuvering case. The command U; is a speed command in the cross range 


direction and not a bearing rate command. While the former is a constant, the 
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latter is range-dependent. Similarly the range rate command U, reflects an as- 
sumed constant range rate for the nonmaneuvering case. This again is an incorrect 
assumption for a target moving along any straight line other than the line of sight. 


The impact of these assumptions will become more apparent in the next section. 


C. MULTI-MODEL ESTIMATION 
1. Concept 

Use of the classic Kalman filter as an observer for a maneuvering target 
is possible only if provision is made for the unknown maneuver command. The 
method suggested by Sin,;e. does not suffice in itself since it dictates maintaining 
high Kalman gain; this inhibits effective measurement noise filtration. 

A better procedure is known as the Multi-Model (MM) [43] and was first 
applied to the MP tracking problem by Moose. In the MM a number (N) of 


commands are hypothesized to form a command bank vector UI defined as: 
Ur= (U,, U2 enOy ). (2.14a) 


The bank is ideally centered around the mean command value and the commands 
in the bank are evenly spread to span the full range of possible commands. A 
corresponding bank of Kalman filters (models) is formed and each filter in the 
bank separately tracks the target using its particular hypothesized command and 
the common measurement Z. The estimated states X; from all the filters in the 
bank are combined to form a conditional mean which is used as the overall estima- 
tor output. The conditional mean depends on the discrete command probability 


distribution represented by the weight vector W defined as 


W =(W, W, ...Wy)? (2.14b) 
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where the weight W; is the conditional probability (weight) that a command U; was 
exercised in the previous time period given all past observations z,. The expression 
for the estimated state X is es X, W;. 

The discrete distribution of W can be computed recursively if the com- 
mand is modeled as a semi-markov process with a known probability transition 
matrix 9. The conditional probability of each command, given the past position 
measurements and the hypothesis, is computable using the Kalman filter’s error 
propagation matrix P. A diagonal matrix A with elements a[z,2] set proportional 
to the conditional probability of the innovation given the hypothesis 2 is also used 
in the recursion which is given below. An overall block diagram of the resulting 
estimator is shown on Fig. 2.2. The development of the basic estimator is reviewed 
in detail in Ref. 24 and will not be repeated here except for the resulting estimator 
equations. Emphasis in the following sections will be directed towards the more 
recent evaluation and modification of the algorithm and to its application to the 
IH multipath tracking problem. 

2. The Multi-Model Estimator Equations 

The equations are written in a generic form for a 3 x 1 state vector X 
representing either range R or bearing B. The specific state vectors, observation 
scalars and the corresponding matrices can be substituted to develop either the 
range or bearing equations. The depth is estimated by means of a first order 
smoother (similar to Eq. (2.10a). The variables and their definitions are given in 


Table 2.2. 
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TABLE 2.2 
MM VARIABLE DEFINITIONS 


3 x 3 state transition matrix. 

3 Xx 1 control gain vector. 

3 X 1 Singer process noise gain vector. 
3 x 1 observation matrix H.= [1, 0,0]. 
3 x 1 Kalman gain matrix. 


Zz 
3 x 1 state vector X = E | at time n. 


3 xX 1 estimate of filter 2 in the bank at time n, 

given past observations 2), 22... 2Zn2 

The observation scalar. 

Number of models in the bank. 

1 x N a vector of 1’s. 

3 x N state bank mairtx. XI[z 7] is the i- th component 

of the state vector Xj of the j-th model, 

i.e., XI = [Ky |X| --|X~]. The matrix is time 

index like X; above. 

The assumed command of the 7*" filter in the bank. 

1 x N command bank, evenly spanning the Umin — Umaz 

interval of allowable commands, UI = [U, U2,...Un]. 

Minimum command in the bank. 

Maximum command in the bank. 

Command bank separation Umaz — Umin. 

1 x 1 command bank step size UI[2] — UI[1]. 

Center of command bank (Umaz — Umin )/2. 

1 x 1 zero mean gaussian measurement noise v ~ N(0, 02). 

1 x 1 measurement noise standard deviation. 

1 x 1 Singer white process noise input. 

1 x 1 standard deviation of the conditional innovation 

distribution given the hypothesis. 

Standard deviation of the command quantization error. 

3 x 3 command quantization error matrix. 

N x 1 hypothesis command probability (weight) vector. 
= (W, W, ... Wy] 

N x N command Markov probability transition matrix, 6 

- the probability that command 7 will change to manele Zin one step. 

dig = ties P.) - os a: where P, 1s the probability of 

unchanging command. 

N x N conditional innovation probability matrix. 

Smoothed output of the estimates X, U,W 

Smoothing coefficient 


ol 


The estimation equations are presented in the order of their use in the 


recursion. The Kalman filter equations are 


XIyin-1 = @* XI p-apn-1 +P: UI (2am 
Paint = $Pa-1jn-16 +Dut+ YoU W (2.16) 
Ga = Pryger YAP een tial (2.17) 
Prin = (I — GnH] Prjn-1 (2.18) 
The. = XIpjn-1 + Gn (en: N-H- XIyjn-1) (2.19) 


Note that the state of the multiple-model XI is the 3 x N state bank matrix whose 
columns represent the states of the individual filters in the bank. The corresponding 
command UTI is an N x 1 row vector which represents the entire command bank. 


D,, is the command quantization error matrix and 1s computed once as 


D, = = i (2.20) 


Since Eq.(2.15) through (2.19) are common to all the filters, the classical single 
channel Kalman gain vector G,, can be computed once and used for all the filters. 


The adaptive command estimation equations are 


jy = HP ini? + 0% — 
—(zn-N—-H-XI[1,1])? /207,, 


€ i =7 
Alijj, = rf uy (2.22) 
=(NA@W,_1]° AOW,-1 (2.23a) 


Wa = dw Wr-1 + (1 — aw)W, (2.235) 


O2 


where the vector of conditional command probability W is recursively computed. 
The conditional command probability W is then used to compute the conditional 


mean estimates X and U of the state and command banks XI and UI respectively. 


a 


U, = UI-W,, (2.24) 


X, = XI,-W,, (2.25) 


Additional first order smoothing is applied to the MM estimate in the form of 


(Xop)n = Ge(Xop)n-1 + (1 —az)Xn (2.26) 
(Uop)n = Gu(UVop)n—1 + (1 — au )Un (2.27) 
(Wop)n = Gu(Wop)n—1 +(1— au) Wa (2.28) 


where the subscript op stands for “output”. 
A high initial Kalman gain is insured by setting the initial conditions as 
follows: 


Tol) ere bf 








XIo 10 (1, 2 = 20 (2.29) 
1 
XIojo (2, 2| => 7 (24 = 29) (2.30) 
XIoj0[3, ] = 0 (2.31) 
Poo(1, 1] = 03 x 10 (2.32) 
o2 
Pojol2,2] = 22 (2.33) 
10 x o? 
Pojol1, 2] = Pojo[2, 1] = ——® (2.34) 
1007 
Porto ee 3] = Pojo 135 1] —_ 72 (2.35) 
Pojo[3, 3] = 0 (2.36) 
1 
Wo = Wil i aes (2:37) 
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3. Estimation Factors 
Some factors in the estimation process emerged as dominant when the 
observer performance was evaluated over a variety of conditions. These factors are 
reviewed in this section. 
a. Hypothesized Commands 
The use of a conditional mean between the set of hypothesized com- 
mand parameters in the estimator leads to the dependence of the estimate on 
the validity and accuracy of the hypotheses and the nature (distribution) of the 
measurement noise. This effect can be demonstrated by the following simplified 
example. Consider an unknown scalar parameter U which is to be estimated from 
a single observation z = U +n where n is a random variable with zero mean normal 
distribution n ~ (0,07). If one assumes a discrete model for U with two equally 
likely possible hypotheses U1 and U2 and U2 > U1, then the Bayes mean-square 


estimate is given by the conditional mean 


j= Ble} = Us Pils) = You, FD POD) (2.38) 


=] =1 
where f is the probability density function of z. This estimate is shown in 
appendix B to yield 


. U, Ole 
U=U,+ > tanh Pex on U.) (2.39) 


where U, and U, are the center and separation of the hypothesis command bank 


defined by 


_ Uy +0, 
= 


Copa) « (2.41) 


U, (2.40) 


o4 
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Fig. 2.3. Conditional mean estimate. 


The dependence of the estimates on U, and U, and on the noise distribution vari- 
ance o is clearly seen in Eq. (2.39). The function U(z) defined by Eq. (2.39) is 
plotted in Fig. 2.3. The bias E = £ {a —U \ of the estimator is given by: 


U, 
207 





é= B10. + Pant | 


; (Ut+n— v-) 2 u | (2.42) 


which after rearranging and extracting the constants from the expectation gives 





Us U, 
E=U.,-U+t+ 3 & {tanh = (UV+n— u.) ; (2.43) 


Given that U. = U the conditional bias can be shown to be zero as follows. Sub- 


stituting U = U, into Eq. (2.43) gives the bias as 


E= SE {tank ( = n) : (2.44) 





OO 


Since the tanh(z) function is odd and the density function of n is even, the expec- 
tation in Eq. (2.43) is zero. 

However the bias is non-zero when U, # U, 1.e., when the hypothesis 
bank is not centered around the actual value. This can be realized for example by 
investigating Fig. 2.3 for values of U larger then U2 or smaller than U,. If U;. #U 
then U will be incorrectly estimated with the amount of error depending on the 

deviation U — U,. 

. This simple analysis which can be extended to the MM case, provides 
a new insight into the operation of the MM estimator and its bias. If the command 
bank is to cover all possible commands of a target which could move at full speed 
both incoming (decreasing range) and outgoing (increasing range), its center should 
be set to zero speed. While a choice of U, = 0 will produce good results when 
averaged over all possible target scenarios, it will produce a bias E for every specific 
scenario with average command E{U} #4 U,. As long as the target is moving in 
one direction (inward or outward) the center of the hypothesis bank will be very 
different from the average command. This will give rise to an estimation bias. 

In the example shown in Fig. 2.4 the command U of a target driven by 
a true command of 10 m/sec is estimated using a bank of 7 filters evenly spanning 
the range —15 to +15 m/sec centered around 0 m/sec. The resulting steady state 
command bias is 0.8 m/sec as shown in Fig. 2.4a. The position error resulting 


from the command bias is 150 m (Fig. 2.4b). 
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Fig. 2.4. Command estimation bias. 
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b. Basic State Separation 
The dependence o. the MM performance on the hypotheses separa- 
tion warrants some further investigation. The measurement and the evaluation 
of the conditional mean are performed in the position z(z) domain. The system 
dynamics as represented in the transition matrix and the Kalman filter translate a 


hypothesized command to a hypothesized position. In particular, consider a system 


X, = @X,-1 + TU + Vu,-1 (2.45) 
with an observation 
Zn, = BX + vy, (2.46) 
modeled by 
X, = $Xp-1 + TUy + Vun-1 (2.47) 


Uy represents the hypothesized constant command here assumed to be different 
than the actual constant command U. The observation noise v and process noise w 
are both normally distributed with v ~ N(0, 02) and w ~ N(0,02,). The variables 
@,T and H are the transition, control gain, and observation matrices respectively. 


The goal is to find the mean steady state deviation & of X from X, i.e., 
€= lim E {Xnin ” Kon} (2.48) 


resulting from the mismatch between the model (Uj) and the actual system (U). 


The Kalman prediction 1s 
Bes > Xi n-1n=1 te LU (2.49) 


and its average is given by 


E{Xnin-t} =$E{Xn-ajn—1} +TUn. (2.50) 
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The Kalman estimate Xyjn 18 
Xajn = (1— GnH) Xajn—1 + Gn2n 


which after substitution of Eq. (2.45) and (2.46) becomes 


Xain = (1— GnaH) Xpin-1 + GnH ($Xn-1 + TU + Wwn_-1) + Gnrn. 


Taking the expectation of Eq. (2.52) yields 
E{Xnin} =(1-G,H)E {is \ + GH¢E {X,-1} + G,HIU 
+ GHVE {wy-1} + GnrE{un} 


and since E{wn} = E{vun} =0 Eq. (2.53) becomes 


E{Xnin} =(1- GaH) E{Xnin-1} + GpHGE {Xn—1} + GaHTU, 


If Eq. (2.50) is substituted into Eq. (2.54) the result is: 


(2.51) 


(2.52) 


(2.53) 


(2.54) 


E{Xnin} = (1- GaH) ($E {Xn—1jn—1} +1Ux) + GaH ($E {Xn-1} + TU). 


Now taking the expectation of both sides of Eq. (2.45) gives 


E{Xn} = ¢E{Xn-1} +TU + VE {wn-1} 


and since E{w,} = 0, 


EB De = oE {Xn-1} soe 


(2.55) 


(2.56) 


(2.57) 


Subtracting Eq. (2.57) from Eq. (2.55), rearranging terms and substituting €,, 


defined as 
En = E{Xnin -Xn} = E{Xnin} — E{Xn} 


o9 


(2.58) 


results in 
E {En} a (I a G,H) ~E Serena + (I al G,H)l (Uy a U) : (2.59) 


At the steady state we can substitute € = €,, for both €, and €,-1 


as well as G = Gq for Gn, in Eq. (2.59) and therefore obtain 

€ = (I- GH) ¢€ + (I—- GH) (Uy - VU) (2.60) 
or alternatively 

(I — (I —- GH) g|&€ = (I- GH) T(Ux — UV). (2.61) 


The desired transformation of the speed command hypothesis deviation Uy — U 


to the steady state average state deviation is thus found to be 
€ =(I-(1- GH)d@|i’ (I-GH)T W@W, - FD) (2.62) 


which has the simple form 
€ = F(Uy —U) (2.63) 


where 


F = (I-(I- GH) ¢]"* (I- GH). (2.64) 


In the MM the hypothesized commands in the bank deviate from the 
actual command value. This command deviation translates into the deviation of the 
corresponding states from the actual state according to Eq. (2.63). The separation 
of the states in the bank is therefore dependent on the command separation, the 


system matrices and the steady state Kalman gain. 
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The dominance of the Kalman gain in Eq. (2.64) is clearly evident. A 
simple scalar example will demonstrate the effect. Consider an attempt to estimate 


a slowly varying constant x defined by the model 
In = a2n-1 +U + wWn-} (2.65) 


where w is a process noise w ~ N (0, 02,) and a has a positive value close to but 


less than 1. The observation is defined by a 


Zn = Hon tn (2.66) 
where H = 1 and v is the observation noise v ~ N (0,02). If the model 
Taln—1 = GZn—1|n-1 + Uy + Wn-1 (2.67) 


is used for the Kalman filter estimate then, the model mismatch Uy —U will give 
rise to a steady state position deviation. Substituting ¢ = a,’ = 1,H = 1 into Eq. 


(2.62) gives the steady state position deviation 


E 1-9 - (Un -U) (2.68) 


~ 1=(1=9) 
where g is the positive scalar steady state Kalman gain which depends on the vari- 
ances ratio of the process and the measurement noise. The deviation, as indicated 
by Eq. (2.68), will vary from a minimum of 0 (for the highest gain case of g = 1) 
to a maximum of ae ( for the lowest gain case of g = 0). 

For the MM the ¢,I, H matrices and the allowable range of maneuver 
commands are assumed known. Under such conditions the only control the MM 
designer has over the state separation is the Kalman gain. A high gain gives heavy 


weight to the new measurement which is common to all the filters; this tends to 
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keep the filters states close together since their differences are not emphasized. A 
low gain gives higher weight to the different individual channel predictions thus 
allowing the states to grow apart. The states diverge until their average innovation 
is large enough to compensate via the Kalman gain for the command mismatch. 
An example is given in Fig. 2.5a showing the positions of the two 
extremes and the center filters in the bank. The Kalman gain (position component) 
was G(1) = 0.073 when the position deviation of the first filter was 1.4km. Fig. 
2.5b shows a similar case with higher gain G(1) = 0.041. Note that the separation 


between the first and center tracks is increased to 2.9Km (note the different scale). 
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Fig. 2.5. State position separation. 
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In Fig. 2.6 a 3-D plot showing the weight vector W as a function of 
time. A sharp and well defined peak (a) leads to a less noisy and overall better 
maneuver tracking performance. This is the case if most of the measurements z 
(Eq. (2.13)) fall within the interval spanned by the state bank XI, and if the 
innovation variance oy is of the order of the position separation. If many of the 
measurements fall outside the interval spanned by XI, as is the case when the gain 
is too high, then a less well defined peak is formed. Such a case is shown in Fig. 
2.6b where the resulting “hesitation” and noisy nature of the weight vector may 
eventually lead to loss of track. 

On the basis of the forgoing discussion, it can be seen that in order to 
minimize the estimation bias one should take the following actions. First one has 
to ensure that the state bank spans the expected spread of position measurements. 
This can be done by tuning the Kalman gain. Secondly one should maintain the 
center of the hypothesis bank close to the actual command. This can be done by 


adaptively relocating the command bank center U, defined as: 


Os 


alr 
Me 


U; (2.69) 
1 


: 
around the estimated command value. Adaptive recentering has not been at- 
tempted in previous work with the MM tracking filter. 

Trimming the Kalman gain, is done by means of an optimized correc- 
tion factor D,. Athans and Chang [4] refer to the optimization as “more of an art 


then a science.” and Ref. 9 attempts to derive an analytical expression for D,,*. 


* The argument made there is that each model in the bank is subject to an 
additional command uncertainty in the interval +40 around the discrete value. 
This is of limited validity since the command uncertainty for every channel is over 


the entire interval Umaz — Umin. However the Singer noise component which is 
added to the D, compensates by maintaining a high enough Kalman gain. 
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Fig. 2.6. Weight vector history. 


Gain control by means of the D, factor reduces the need for the Singer state com- 


ponent w’, which was devised partially for the same purpose. Our approach, which 


was found to be more effective is discussed in Section D. Adaptively recentering 


the command bank, was a more involved procedure than it seem to be at first. The 


method is also described in Section D. 
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To complete the discussion Eq. (2.63) is now extended to case of a 


complete MM state bank with N channels. Recall from Table 2.2 that the state 


bank matrix has the form 


XI = (Ki |Xe]...|Xw] 


(2.69) 


where X,X2...Xwn are the state vectors for the individual filters, and that the 


command bank vector has the form 


UI = [(U;,U2,...Uy] 


(2.70) 


where U,,U2,...Uyn are the corresponding commands. The state bank deviation 


matrix defined as 
EI = (€1, €2,...En] = (Ki — X, Xo — X,... Xn — X] 


is given by 
EI=XI-X-N 


where 


N =[1,1,...1). 


It follows from Eq. (2.63) that the error for the 2*” filter is given by 
€& = F(U; —U) 
and thus €I is also given by 


EI = (F(U,; —U),F(U2 —U),...,F(Un —U)] 


€Il = F(UI-U-N) 
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(2.71) 


(2.72) 


(2.73) 


(2.74) 


(2.75) 


(2.75) 


namely 


XI —X-N=F(UI- UN). (2.76) 


Now if each control U; is changed by the same amount AU, then the change in the 
entire command bank AUI is AU -N and the new command vector will be given 
by UI + AU-N. From Eq. (2.63) the corresponding change in the state bank 
matrix is given by 


AXI = F-AUI=F-AU-N. (2.77) 


The foregoing discussion implies that for an evenly spaced command the state 
bank matrix is also evenly spaced along the positior -:-eed and acceleration axis, 


the state spacing AX is given by 
AX =F.- AU. (2.78) 


Eq. (2.78) compared very favorably with simulation results (the predicted deviation 
was within +0.01% of the simulation result.) and turned out to be very useful in 
ridding the MM tracker from one of its inherent biases ( in Section D). 
c. Smoother Order 

A first order averager, Eq. (2.26), (2.27), and (2.28) was added to 
smooth the filter output in order to reduce the output variance. The price paid for 
the smoothing is a tracking lag error for targets of constant speed [9]. This results 
from the mismatch between the target model (second order) and the averager (first 


order). The steady state error E can be derived using the final value theorem 
. z£—l 
len = lim ae) (2.79) 


where X(z) is the z transform of z,. Consider the scalar, constant rate target 
position Z, given by 


Cn = Uy,- P(N). (2.80) 
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_ where u(n) is the unit step function and U, is the rate parameter [position units/sample 


time JT], and a smoothing filter of the form 
Yn = 2Yn-1 AF (1 — NL . Q2.81) 


The observation z is here assumed to be noise free so that z, = tn. The 


smoothing filter produces a lag error E [position units] defined as 


a — Un — er. (2282) 


Sa 


whose value at steady state is sought. The analysis in the Z domain proceeds as 
follows. Let H,,) be the system response and 4(,),,) be the Z transformations 


of z, and y,. Then 








d Upoz 
Me) = —2- Fo Uw] = Cap (2.83) 
1—az! 
= ———_—__ 2.84 
Ui) i ( 8 ) 
Uja 1 
Ez) = Vz) — X2) = (Ae) — IX) = + T- (2.85) 
Applying Eq. (2.79) yields the steady state error | 
; z—1 aU aU : 1 
ao z Gaile Decck—oe 2 2 age. 
and the steady state lag error is thus 
a 
l-a 


indicating large errors for heavy smoothing (a close to 1). 

In the time domain Eq. (2.81) averages the past estimated value 
X(n —1) with the current measurement of z(n). This produces an error even 
for the case of no measurement noise since the position z is different at the two 
successive time samples. This intuitive understanding of the source of the error led 


to the development of the improved advancing smoother discussed in Section D. 
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d. Coordinate Linearization and Decoupling 

An important aspect of the coordinate decoupling and linearization 1s 
that the target helmsman commands in the decoupled system are attached to a 
moving coordinate system, namely the target range and cross range. The drawback 
of such a system is that a (nonmaneuvering) target moving with constant velocity 
along any straight line other than the line of sight has to be described as maneu- 
vering since its range and cross-range speeds are constantly changing. However in 
the MM the probability that the command will not change is assumed to be very 
high (P, close to 1). 

These two opposing assumptions produce an error which is more pro- 
nounced at short ranges where range and cross range accelerations are higher. At 
short ranges the normalized range rate (range rate/range) is also large. This im- 
plies that there is a need to update the bearing channel transition matrix more 
frequently. Under the foregoing condition the advantages of this linearization and 
decoupling approach are marginal. 

A potential remedy to the coordinate linearization problem may be 
to return to an XY Cartesian system or to adaptively modify the probability 
transition matrix @ to reflect a predicted linear motion. These ideas were not 
pursued further however in order to concentrate effort on the IH effects of the 
MP tracking. 

e. Measurement Noise 

In passive localization, TDOA from the target via different acoustical 
paths to separate parts of a sonar array are converted to depth, range and bearing 
information. The process is error prone both in the time delay estimation and in 
the transformation to target position (range depth bearing). This is one of the 


major topics addressed in this research and most of Chapter Three is devoted to 
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this problem. However the noise dependence on range affects the evaluation of the 
target tracker itself, and so is described briefly here. 
The time delay estimation error depends on the signal to noise ratio 
(SNR) of the received signal. The SNR in turn depends on the transmission loss 
and therefore indirectly on the range. Thus, the delay noise is range dependent.* 
In addition, the time delay estimation limits the estimated delay to 
positive values. An improved noise model is used for the tracker performance 


evaluation. This noise model is described in the Section D. 


D. IMPROVEMENTS AND MODIFICATION ” 
The following modifications of the tracker were introduced to improve it and 


evaluate its performance with respect to the issues discussed in Section C. 
e An adaptive instead of a fixed hypothesis bank was incorporated. 


e A second instead of a third order system was incorporated in the model. The 
second order system was demonstrated as sufficient for the target when the 


steady state Kalman gain is adjusted by the correction factor D, (see eq. 
(2.16), (2.17), (2.20)). 


e A new second order smoothing algorithm was developed and used to replace 
the first order smoother. 


The modified tracker was than evaluated with 
e A more realistic range dependent and nonnegative noise. 


e Emphasis on range bearing coordinate decoupling and linearization at short 
ranges. 


* With this limitation in mind the disadvantage of the coordinate linearization, 
discussed above, at short ranges is further emphasized. 
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1. Adaptive Command Bank 


a. Concept 

The problem addressed by this modification is the bias error arising 
from a hypotheses bank which is not symmetric around the actual command. A 
possible solution to this problem, cited earlier, is to adaptively recenter the hypoth- 
esis bank around the actual value. Recall that when the center U, of the command 
bank (referred to here as “the center” for convenience) is equal to the value of 
the true command, this bias is removed. Since the actual value of the command 
is not known, our algorithm uses the command estimate U,, instead. The main 
idea is therefore to periodically recenter the command bank around the estimated 
command value. 

A straightforward implementation is to average the estimated control 
over a time interval long compared to the system and the observer time constants 
and to feed back the average as the new center. Specifically the difference between 
the estimated command Ug, and the center U, is averaged using an autoregressive 


filter to produce the center deviation U¢q. 
Ved noe aU ed n—1 + (1 7 a)( Vag n—l] — Uy =) (2.88) 


This average value is used to dynamically shift the center and the complete com- 


mand bank. 


oe n+1 = Coe tin (2.89) 


Wilkes: |: =n ietialad eee (2.90) 


Although this method was implemented, its performance was poor for the following 


reason. The combination of the command hypothesis UI and the weight vector 
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W is a discrete description of the command probability distribution estimated by 
the Bayesian recursion (Eq. (2.23)). Updating the command hypothesis vector but 
keeping the weight vector unchanged effectively modifies the command distribution 
and disrupts the recursive command estimation. 

Further, the command bank yields, over some time, a corresponding 
bank of hypothesized states in the form of the matrix XI which based on Eq. (2.76) 
is 

XI = |€,&...En| + X-N=F(UI+U -N). 
Different states correspond t« the new updated command bank. Shifting the com- 
mand hypotheses bank introduces an undesirable transient in the state vectors of 
the filters. This transient seriously degrades the tracking. Furthermore, the feed- 
back of the state vector’s transient into the recursive command estimation process 
(Eq. (2.22) and (2.23a)), may lead to instability of the entire tracker. 

It is clear therefore that updating the command bank center alone 
disrupts the recursive estimation process unless there is a corresponding update 
of the weight vector W and the state matrix XI. In simulation recentering the 
command bank without updating the weight vector and the state matrix resulted 
in complete divergence of the command estimate. 

b. Complete Model Updating 
An alternative approach was devised where the complete model was 


updated and the instability was eliminated. The method was as follows: 


e The update was restricted to occur at discrete events in time when the averaged 
deviation of the center U, from the estimate U,, reached a preset threshold. 


e The commands and the command center were restricted to take on a set of 
discrete quantized equally-spaced values. Updating was done by shifting the 
entire bank up or down one position. 
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e The complete model including command and state banks, and the weight 
vector were updated. 


The time discretization was intended to prevent the continuous dy- 
namic feedback which led to instability. The time constant used for the exponential 
command deviation averaging was set larger than both the system and the MM 
estimator time constants. This was done to ensure complete transient recovery 
from one update before another update is allowed. 

The model update was implemented in a cyclic manner. The least 
likely filter in the bank, i.e., the one corresponding to the command furthest away 
from the estimated command, is dropped. A new and more likely command is 
added at the closer end of the bank. The probability of the new channel is initially 
set to that of the one that is dropped, so that the total probability 3 W; 1s kept 


=] 


equal to 1. If the update takes place at time k then 
UI,+ = UI, t AU-N (2.91) 


where the subscripts k~ and kt represent the control prior to and after the update 
and the sign of AU depends on the sign of the averaged center deviation U.q (Eq. 
(2.88)). 

A block diagram of the modified MM is shown in Fig. 2.7 and Fig. 
2.8 shows an example of a model update process assuming the update takes place 
at time index n = k. The variables at time k~, i.e., just prior to the update, 
are shown in Fig. 2.8a; the variables immediately after the update (time kT) are 
shown in Fig. 2.8b. Prior to the update the bank is centered around U, = 0 
and spans the interval —15 to +15 [m/sec] as shown along the horizontal axis. 
The second (lower) horizontal axis represents the first component (position) of 


the state vectors of all filters in the bank (this is the first row of XI) before the 


2: 


update. This spans the interval 8.5 to 11.5 km. During the update the first 
hypothesis (U;),- =~—15 m/sec is dropped and a new hypothesis (U7 ),+ is added 
at +20 m/sec with a corresponding state position XI{1, 7],+ at 12 km (Fig. 2.8b). 
The weights W; (conditional probabilities) for all the filters except that of the new 
one are shifted one position. This is done in order to maintain the correspondence 
with the hypotheses. The new filter assumes the weight of the filter that is dropped 
(0.01 in the example shown). The command bank is now properly centered around 
5 m/sec. The command estimation transient is shown in Appendix C to be given 
by 

Vie en = NAUI(VADE- (2.92) 


For the specific case in Fig. 2.8 this transient is therefore 
Utransient = 7-5- 0.01 = 0.35m/sec. . 


This minor transient is further smoothed by the output smoother and its effect on 
the output Uo, is negligible. 

The command in the bank and the center U, are restricted to the 
evenly spaced discrete values. The center update is further limited to a one position 
shift up or down from its current position on the discretized scale. The command 
bank update is. thus implemented by means of the linear shift given by Eq. (2.91). 
Since the states of the bank are also evenly spread over the state scales (position, 
speed and acceleration). The corresponding update of the states matrix XI can 


thus be implemented by means of a linear shift as follows 
XI,p4 = XI,_ + AX -N (2.93) 


where AX is the constant vector F- AU - N (Eq. (2.78)). The complete model 1s 


thus updated by a simple one position shift. The shift is linear for the command 
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Fig. 2.7. Model update block diagram. 
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UI and the state XI, and circular for the weight vector W. The simplicity of 
the model update procedure (which requires no computation) resulted from the 
constraint of the center U. to assume values on the discrete command scale. 
Tracking with the complete adaptive model (updating UI, W and XT) 
is shown in Fig. 2.9. Notice the command and position biases developing just prior 
to the update 2.9c which is reduced by the two updates at the 18'* and a5 
minute, and the almost unnoticeable transient. Further persistence of the same 
target motion would lead to additional updates which would eventually reduce the 
bias completely. 
2. Second Order Target Model 
In Eq. (2.8) the second order target model was augmented to a third order 
model in order to account for the uncertainty in the maneuver command. Prac- 
tically this augmentation effected the tracking by increasing the error covariance 
matrix P which in turn increases the Kalman gain (Eq. (2.14)). Since the Kalman 
gain is so critical in the MM estimator it is additionally trimmed by the factor Dy. 
In the original third order Singer model the third state variable w was 
the only state variable that accounted for the target acceleration. In the MM 
the command (and the corresponding acceleration) is separately estimated as the 
independent state U. With both the acceleration estimate and the Kalman gain 
handled otherwise the use of the third state w is redundant and can be eliminated. 


Elimination of this state results in the following second order model 
XI r= $, XI, n=l PoUR-1 (2.94) 


Where XI,H 2 and Il are the first two rows of XI,H andT respectively, and ¢, 
is the 2 x 2 upper left block part of ¢ (see Eq. (2.15)). This configuration will save 


computation time and simplify the gain adjustment procedure. 


16 





AL: CONTROL * SPEED +; ESTIMATED CONTROL VY CENTER A 


CONTROL [M/S] 


2 





TIE [MIN] 


-—200 —100 0 
COT: a a ee 


RANGE [M] 
-300 


-400 


20 


ee ie St 





RANGE ERROR 


40 
TIME [MIN] 


60 


ADAPTIVE WEIGHTS 
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Fig. 2.10 compares the tracking of a third order model with steady state 
Kalman position gain G[1] = 0.041 to the tracking of a second order model with 
practically the same position gain (G[{1] = 0.042). The gain of the second order 
model is controlled by 7, which affects the gain via Dy, = To2P’. The performance 
of the two models is similar except the third order system has a larger overshoot. 

3. Second Order Smoother 

An improved second order smoother referred to as an advancing smoother 
was introduced [25] in order to reduce the constant speed target lag error (Eq. 
(2.75)). The basic idea is to average the new filter output Z,), with the predicted 
point Zpjn—1 rather than with the previous estimated position Zp-1|n-1. The 


current estimated speed is used to generate the predicted position according to 
(Xop)n = ie ag 0] ; (Xop)n-1 a (1 7 a)(1, 0; 0)X{1Jn. (2.95) 


If carried out to the full extent, this principle would turn the smoother into a 
Kalman filter by itself. The need for this additional smoothing, which is usually 
handled by means of a reduced Kalman gain, results from the optimizing the gains 
to meet the overall MM requirements. This forces the additional external smooth- 
ing. 
4. Improved Measurement Noise Model 

A range dependent, nonnegative noise is used in this simulation in place 
of the range independent normally distributed noise used in past work. 

a. Dependence on SNR 

The dependence of the delay estimation noise standard deviation (STD), 

on the range, via its dependence on the SNR, was modeled here in two ways. In 


the first model the range dependence is lumped into the following equation. 


STD;(R) = ST Dip - RP (2.96) 
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where STD,, is the i-th TDOA error initial standard deviation at 500 m from the 
source, and p is the delay noise range power factor which depends on propagation 
loss. The variable p takes on values in the range of 1 to 3. A low value for p (p = 1) 
corresponds to a high SNR and cylindrical propagation while a higher value (p = 2) 
corresponds to low SNR and spherical propagation. Higher p values corresponding 
to increased loss due to ray bending were not considered. 


In the second model, a more recent and detailed description was used 


Pd 
Theme’ 


based on Ref. 14. Here the exact dependence of the variance on the bandwidth, 
the observation time T and the SNR was implemented. The parameters of the first 
simplified noise model (ST'Dp and p) were set to provide an overall similar noise 
variance which is similar to the results of the second noise model. 
b. Nonnegative Time Delay 

Another feature of the time delay estimation for a single receiver is 
the inability to distinguish a leading signal replica from a lagging one. This results 
from the symmetry of the autocorrelation function (ACF) of real signals. The 
nonnegative time delay effect was simulated by applying the absolute value operator 


to the time delay after the noise was added. That is 
Tom = |%i +; (2.97) 


where v; ~ N(0,STD?(R)). This operation produces a bias primarily for 
sTDi(R) <1, which is also discussed in Appendix D. 

After the addition of the noise the simulated TDOA undergoes a non- 
linear inversion to depth and range which also produces a bias. The bias was 
studied for the idealized, homogeneous and space invariant delay noise, by Moose 
in Ref. 11. Investigation of this error and its correction for the JH case, and the 


means to account for it are discussed in Chapter Four. 
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c. Depth Range and Bearing Noise 


When the filter was evaluated by itself (without the prefilter) noise 


was added to the depth and range measurement directly. Normally distributed 


range dependent noise was used such that 


Dm = Dt+ug 


Rn = R+0, 
where 


vg ~ N(0,STD?(R)) 


ve ~ N(0,STD3(R)) 
and the standard deviations given by 


STD4(R) = STD4o - R?® 


STD,(R) = STDyo - RP 


(2.90) 


(2.99) 


(2.100(a)) 


(2.101(b)) 


(2.101(a)) 


(2.101(b)) 


with STDa,. STD, being the initial STD at 500 m and p, the range dependence 


factor. 


Range dependent noise was added to the bearing channel as well in 


order to maintain a consistent uniform simulation. The bearing measurement is 


thus 


B, =Bt+uwv, 
where 


v, ~ N(0, STD?). 
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(2.102) 


(2.103) 


The bearing error standard deviation used here is given by: 
SOC sg — ope it. (2.104) 
where ST Dj, is the bearing error standard deviation at 500 m from the source. 


E. SIMULATION 
1. Model Description 
The simulation model was designed to evaluate the overall MP tracker 
performance, both in the idealized homogeneous case and later in the realistic and 


IH environment | Unapter Four). The following were set as subgoals for the design: 


e Support of a multilevel system evaluation starting from the individual algo- 
rithm up to the complete integrated system. 


e Provide a detailed evaluation of the 3-D maneuvering target tracker perfor- 
mance with emphasis on the 3 axis integration and the tracker modifications. 


e Support of a gradual departure from the idealized homogeneous MP config- 
uration towards the more realistic IH case. The transition includes both the 
simulated IH measurement and the algorithm designed to deal with this dis- 
tortion. This again is in support of Chapter Three. 


In order to maintain a consistent and unified simulation environment the 
model was designed to meet all of the above subgoals by optional use of its various 
components. The model is discussed here in its entirety; however the details of 
some of its components are addressed in the following chapters. 

A block diagram of the simulation model is shown in Fig. 2.11. The 
simulation proceeds as follows. Depth, range and bearing target position D, R, B 
is generated by the scenario generator. The. depth and range are converted to 
TDOA 7,72 by the MEDIUM. Noise generated by the NOISE module is added to 
the TDOA to form the measured 11m72m Which are fed to the PREFILTER. The 
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prefilter performs a functional inversion to compute the depth and range measure- 
ments Dm Rm. Noise is also added to the bearing to simulate a realistic bearing 
measurement B,,. The TARGET TRACKER uses the three noisy measurements 
to produce a 3-D target estimate. The estimate includes position R,p, Bop, rate 
Roy, Bop and command Uy op, Us op, for the range and cross range axis, and posi- 
tion D,, for the depth axis. The noisy measurements and the tracker estimates 
are compared to the output of the scenario generator to produce the measurement 
and estimation errors. Details of each of the blocks are described below. 

The SCENARIO GENERATOR contains a model for both the target and 
its pilot. It includes the target dynamics (Eq. (2.10)) as well as a prescribed target 
trajectory and it outputs the target position in depth, range, and bearing. Three 


different types of pilots were programmed. 


Pilot 1 describes the target motion by an initial state and a set of discrete 
and independent sequences of maneuvering commands along each of the 
three axis. 


Pilot 2 is similar to Pilot 1 except it provides a continuous gradual change of 
the maneuvering command. 


Pilot 3 prescribes a target moving along a straight horizontal line with constant 
speed (Eq. (2.8)). As discussed in Section C this target will be seen as 
maneuvering by a cylindrical coordinate tracker. 


The MEDIUM converts the target depth and range coordinates to time 
delay differences 71,72 (or t,t2). It includes a model for the effects of the medium, 
the receiver, and the noise. Two groups of models were programmed, namely the 


homogeneous group and the inhomogeneous group described below. 
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Fig. 2.11. Simulation block diagram. 
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The homogeneous group. 


Medium 1 This is a bypass mode where the depth and range are used directly, 
without the conversion to time delays and back to depth and range. 
This mode is used to evaluate the MM tracker by itself. 


Medium 2 The approximate formulas valid for R >> D, [Ref. 9] 


_ _ 2Do(Do - D) 
a RC 

24D, BD, (a, — Dy + D) 
i a = 7 ca 


(2.105a) 
(2.1050) 


were used for the direct function. 


Medium 3 The exact c}: sed form analytical expression of time delays as func- 
tion of depth and range derived by Hassab (Eq. (1.1)) were used here. 
The implied assumptions here, besides straight line propagation, are 
that the delays are perfectly resolvable and can be associated with their 
corresponding paths. 


The inhomogeneous group. 


Medium 4 This is the only medium that was implemented in the group. Here 
the three assumptions implied above were removed and an IH medium 
with finite resolution and non-associable time delays is generated. The 
theory and implementation of Medium-4 are the subject of Chapter 
Three. TDOAs for the realistic IH medium are noted as tj, te. 


The NOISE models were used for the TDOA measurements. They were 


noise type (NT) one and two. NT 1 used Eqs. (2.96) and (2.97). NT 2 used Eq. 


(2.96) and the more elaborate noise model which provides the STD; as a function 


of range and the target signal to noise ratio at range of 1 m (SNRy). The model 


and the parameters it uses are described in Appendix D. 


When no medium and prefilters were used, noise was added directly to 


depth and range using Eq. (2.98) and (2.99). Noise was also added to the bearing 


according to Eq. (2.102). 
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The PREFILTER converts the noisy TDOA estimates transforming them 
to depth and range. This is done by inverting the nonlinear divect function. Six 
different prefilters were designed and grouped into two groups: homogeneous and 
inhomogeneous prefilters. 

Homogeneous prefilters. 


Prefilter - 1 This is the bypass case where no functional inversion is done since 
the inputs are depth and range, corresponding to medium 1. 


Prefilter - 2 This uses the following approximated relation 


Ba (2.106a) 


(2.1068) 


valid for R >> D corresponding to medium-2 [Ref. 9]. 


Prefilter - 3 This is the exact analytic closed form relation of Eq. (1.2) corre- 
sponding to Medium 3. 


Inhomogeneous prefilters. 


Prefilters - 4,5,6 These correspond to the inhomogeneous medium, with finite 
TDOA resolution and non associable delays of medium 4. These again 
are discussed further in Chapter Three. 


The TRACKER uses the measured depth, range and bearing as inputs for 
the 3-D maneuvering MM target tracker. Third and second order target models 
with first and second order smoothers were realized. The command bank was either 
fixed or adaptive as discussed in Section D.1. 

2. Description of the Output Plots 

Three main plots are produced as output for each axis of the cylindrical 

coordinate system; they are the position, the speed and the command plots. When 


a combination of plots is presented the plot symbol markings are synchronized in 
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time on all the plots of a given run to provide a means for cross reference between 
them. 

The values of the simulated Runs were recorded and plotted only every 
fifth sample time to avoid cluttering of the picture. An example of a full sample 
rate plot (with a point plotted for every sample time period) versus a reduced 
sample rate (point plotted every fifth sample time) is shown on Fig. 2.12. 

The position plots (depth, range, bearing, and X,Y as functions of time) 
include the true, measured, and estimated positions. The true position is the 
scenario generator output D,R,B; the measured position is the prefilter depth 
and range outputs D,,Rm and noisy bearing Bm; and the estimated position is 
the tracker output Dop Rop Bop. Position error plots show the difference between 
the measured or estimated positions and the true position. 

The speed plots present range or cross range speed as a function of time. 
Various combinations of the actual and estimated speed and speed command are 
included. The center of the command bank is also shown. Note that when Pilot-3 
was used, constant U, and U, are applied (Eq. (2.8)). The commands U, and U; 
are not used. The model however does estimate U, and U;, which are presented 
on the plots. 

Command weights of the recursive Bayesian estimator are plotted in 3-D 
plots of the weight vector as a function of time. These plots provide a good measure 
of the MM adjustments and tracking quality. 

3. Parameters Selection 

The large number of parameters required for each run dictated a strict and 
well controlled mechanism of parameter selection. This was especially required for 
the realistic [H case where there was a large amount of ocean and acoustic data 


and parameters. The data was divided into the following three categories: target 
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Fig. 2.12. Plot sampling rate. 
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- and scenario; medium and prefilter; MM parameters. The detailed definition of 


the parameters and their values in the simulation is presented in Appendix D. 


F. SIMULATION RESULTS | 

Five straight line trajectories are presented differing mainly in the noise used. 

Run 1 used range dependent depth and range noise without the multipath 
measurement (Mdn =1,Pfn=1). The range dependence p was set at 1 and the 
depth and range error STD at 500 m were set to 1 and 50 respectively. 

Runs 2 - 5 used range dependent TDOA noise based on the improved noise 
model (NT = 2) with SNR and p set according to Table 2.3. 

A detailed list of the other parameter settings 1s included in Appendix E. 

TABLE 2.3 
NOISE PARAMETERS 


Run SNR» [dB] p 


2 o0 1 
3 70 Z 
4 60 Z 
+) 50 Zz 


Results of Run 1 are shown in Fig. 2.13. Range (a) depth (b) and range error 
(c) are shown. The increase of measurement and tracking errors with increasing 
range is seen on all three plots. The very effective noise filtration of the filter is 
demonstrated in Fig. 2.13c after the completion of the initial transient (around 
the 12'* min and on). The range measurement noise STD which was 1000 m at 10 
km was reduced to less than 100 m. Range error reaching 17km develops around 


the CPA which occurs at the 17'* minute at a range of 3 km. The results of this 
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Run set a reference performance for the MM maneuvering target tracker by itself 
| (without the prefilter). 

Results of Run 2 are shown in Fig. 2.14. The performance of the MM 
tracker and the homogeneous medium and prefilter are demonstrated with range 
dependent TDOA noise. 

As can be seen, the cylindrical propagation (p = 1) provides very favorable 
conditions even when the initial SNRo is relatively low (50 dB at 1 m range). Note 
the four model update cycles the MM goes through in response to the cheamial 
range rate. Also note the vanishing range tracking error which results from the 
use of the advancing smoother (25*" minute in Fig. 2.14a). Good bearing tracking 
and the combined horizontal tracking are shown for this run in Fig. 2.15. The XY 
tracking error at the CPA is clearly evident. 

The effect of increasing TDOA noise as a result of decreasing target signal 
strength (Run-3, 4) is shown in Figures 2.16 and 2.17. The failure of the —_— 
in Run 4 past the range of 15Km at the 27th minute of this run is also evident in 
the command estimation plot shown in Fig. 2.17c. 

Very noisy TDOA measurements resulting from low SNRog and spherical 
spreading (Run-5) is shown on Fig. 2.18. Large estimation bias which results from 
the nonnegative effect is clearly seen at ranges larger than 12 km. The depth is 


also underestimated at these ranges. 
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Fig. 2.15. Run-2: Bearing and XY plot. 
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Fig. 2.18. Run-5: High range dependent noise SNRo of 50 dB. 
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SUMMARY OF TRACKER EVALUATION 


The performance of the modified target tracker is summarized below. 
Generally the tracker is well suited for 3-D maneuvering target tracking. 


Measurement noise variance is significantly reduced by the tracker. The range 
measurement STD of 1000 m can be easily handled and reduced to around 50 
m. 


The performance depends strongly on range, especially when range dependent 
MP measurements are used. However up to 15-17 km the filter tracks well 
without readjustments. 


Tracking errors that can reach 2 km in range develop at short ranges for 
nonmaneuvering targets moving along straight lines due to the coordinate 
decoupling and linearization. 


Recentering of the command bank around the average actual command is 
required in order to remove the inherent MM estimation bias. The new model 
update technique devised here effectively accomplishes this task by recentering 
the command bank around the estimated command value. 
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lt. INHOMOGENEOUS REALISTIC MULTIPATH DEPTH 
AND RANGE MEASUREMENTS 


A. INTRODUCTION 
The following idealizing assumptions mentioned in Chapter One and used in 


previous broadband MP tracking studies are removed in this chapter: 
e Homogeneous straight line propagation. 
e No propagation and reflection loss. 


e Infinitely resolvable TDOA. 


e Ability to associate delays with acoustical path. 


The impact of these effects on the tracker is investigated in Section B. Section 
C presents the improved prefilter which provides the functional inversion of time to 
depth and range measurement for the IH realistic medium. Chapter Four analyzes 
the performance of the new inversion algorithm as a component in the overall MP 


tracking system. 


B. REALISTIC MEDIUM AND RECEIVER EFFECTS 
1. Propagation Loss and Reflections 
The most dominant effect the medium has on the propagating acoustic 
wave is the transmission loss. Three contributors to the loss are discussed: the 
spreading, the absorption, and the imperfect reflections from the boundaries. Spe- 
cial attention is given to the effect that the reflection from the non-ideal boundaries 


has on the multipaths. 
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a. Spreading Loss 
The spreading loss, as the name implies, is the result of the spreading 
of the acoustic energy over a growing area, as the distance from the source increases. 
In spherical spreading, the energy spreads evenly over the surface of a sphere. This 
makes the intensity, the power per unit area, inversely proportional to the squared 


range. The relation is 


I(R) = Ip - R~? (3.1) 


where I is the intensity at a distance of 1 m from the source. 
Under other propagation conditions different range dependency re- 
sults. For example, with cylindrical propagation the intensity is inversely propor- 


tional to the range itself, that is 
IR) =Ip-R™ (3.2) 


In deep water, away from the surface or the bottom, one finds spherical 
propagation to be a good assumption for short ranges. In shallow water, where 
MP effects are more dominant, cylindrical propagation is a better representation. 

In the IH case the wavefront is distorted by the variation in the speed 
of sound. Although the resulting spreading loss is not simple, a dependence on 
the p-th power of range can still be useful as a first order approximation. The 
approximation 1s 


I(R) =I): Ro? (3.3) 


where the empirical constant p, which we call the range propagation power, is not 


limited to integers. 
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b. Absorption Loss 
Part of the energy of the propagating acoustic wave is absorbed by 
the water and eventually is turned into heat. The absorption loss depends on the 
salinity of the water, the frequency, and the distance traveled. In short to medium 
ranges (below 15 km), and at sufficiently low frequency (below 10 KHz) this is not 
a dominant phenomenon and will therefore not be considered here. 
c. Reflections From the Boundaries 
The laws of pressure and velocity continuity across the boundary gov- 
ern the reflections from the surface and the bottom [Ref. 23]. A wave propagating 
in medium 1 impinging on a boundary at an angle ( is both reflected back at an 


angle @; and refracted into medium 2 at an angle {2 as shown in Fig. 3.1. 





Fig. 3.1. Reflection from a boundary. 
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The acoustic pressure reflection coefficient defined as 
a (3.4) 
Pi 
where p;, p; are the acoustic pressures of the incident and reflected waves in medium 
1. The reflection coefficient is given by 


_ 42008 B, — 2, cos 42 


~ Zz cos B + Z cos 62 — 


where Z, and Zz are the specific acoustic impedances in media 1 and 2, respec- 
tively. The specific acoustic impedance is defined as the ratio between the complex 
amplitude of the acoustic pressure p and that of the magnitude of the particle 
velocity U. The angle #2 is given by 62 = sin’ (@ sin 6, | where C and C2 are 
the speeds of sound in the two media. 

The specific acoustic impedance for a plane wave can be expressed in 


terms of the speed of sound C and the equilibrium density pp of the medium as 
a Po ° C (3.6) 


which is known as the characteristic impedance. 

Since the characteristic impedance of air is much smaller than that of 
water, the surface boundary has a reflection coefficient close to —1 (see Eq. (3.5)). 
Further, since the characteristic impedance of the bottom is generally larger than 
that of the water, the resultant reflection coefficient of the bottom is positive. 

The reflections from the water air boundary becomes somewhat more 
involved when the layer of air bubbles formed near the surface by the wind and 
the ocean waves is considered. But unlike the air (which is isotropic) the bottom 


may be lossy and anisotropic. In particular, shear waves excited by the impinging 
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acoustic wave travel in the bottom with different characteristics than the refracted 
longitudinal waves. These effects give rise to less predictable and complex acous- 
tic impedance and complex reflection coefficients which strongly depend on the 
bottom type. For a boundary with frequency dependent characteristic impedance 
Z2 =1r2 + Jjz2, the pressure reflection coefficient will vary both with frequency and 
angle of incidence, and will depend on the bottom type and sea state (for surface 
reflection). Exact prediction of the phase of a reflected wave is made difficult by the 
large variation in pressure reflection coefficients. This phenomena has a marked 
influence on MP tracking. The intensity, which is related to the square of the 


pressure, has the reflection coefficient ['; given by 
Ty, =|T|’ (3.7) 


which also depends on the above ocean and acoustic wave parameters. In a very 
simple practical analysis, the energy loss due to the various boundary effects is 
lumped into a single number called the reflection loss. The reflection loss depends 
mostly on the sea state, the bottom type, and the frequency band. Loss values in 
the range of 3 to 30 dB per bounce are commonly found. The sensitivity of the 
reflection phase to frequency and angle of incidence has a marked influence on MP 
measurements. 
2. Inhomogeneous Medium 
a. Speed of Sound Variations 


The speed of sound in the water is given by the formula [23] 
C = 1449.2+ 4.67 —0.055T +0.00029T + (1.34 —0.017)(S — 35) +0.016D (3.8) 


where 7 is temperature in degrees centigrade S is the salinity in parts per thousand, 


and D is the depth in meters. 
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The sound speed is obviously not constant in the vertical water column, 
due to change of the pressure and the temperature with depth. The variation of 
sound velocity with depth, referred to as the sound velocity profile (SVP) is the 
main source of inhomogeneity in the acoustic ocean medium. Waves propagating 
in such a medium are refracted in a complicated manner, and obey the second 


order linear partial differential scalar wave equation [20] 


Volts) — Gags * GA eltst) = X(T) (3.9) 


where y(t, r) is the velocity potential at time ¢ and position r. C(r) is the speed of 
sound at position r, and Xm(t,r) is the source distribution (volume flow rate per 
unit volume). The acoustic particle velocity vector U(¢,r) and pressure p(t, r) are 


given by 


p(t, x) = ~po=o (tr) (3.10) 


Ce HHVoer,: (3.11) 


Solution of the wave equation for the general case is very difficult. 
However, when the speed of sound is a linear function of depth the approximate 
ray acoustic solution yields simple closed-form equations. The detailed develop- 
ment of the ray acoustics approximation is based on the WKB approximation and 
is presented, for example, in Ref. 20. The resulting ray tracing equations are 
presented below. 

b. Ray Tracing 
The parameters and geometry of a typical ray-trace are shown in Fig. 


3.2. Note that the rate of change of sound speed with depth is a constant known 
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Fig. 3.2. Acoustic ray path. 


as the sound speed gradient g, which has units of [sec~'], and the speed of sounu 
C is given by 
C(D)=Cyo+g9D (Sere) 


where Co is the speed of sound at the surface. 
In transition from depth D, where the sound speed is C} to depth D2 
with sound speed of C2, the ray angles 8, and {2 obey Snell’s law 


sin By “¥ C1 


is; ee (3.13) 


Thus a ray emanating from a source at D,,,, at an angle (,, reaches depth D2 


at an angle B2 given by 


B. = sin~' (= sin As) = sin’ Been sn 8 (3.14) 


and at a range R2 given by 





Cy | 
R, = Ry + | - B, — cos B2) (3M) 
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which can be evaluated by solving Eq. (3.14) for 62 and substituting into eq. 
(3.15). The travel time T and the arc length S of the path are given by 


rai tan( {2 /2) 





7 tan(B; /2) © 
and 
Co 
= Pe, (Bo — Bi). KGisticig) 


The ray traverses the depth range (DR) plane along of a circle with radius Ryay 


given by 
Co = C1 COS By 


Reavy = 
ray g g 


(3.18) 


centered at the initial range and the conceptual depth Deenter. This is the point 
where the speed of sound would have become zero if the gradient had been constant 


to that point; it’s value measured from the surface is therefore given by 


Deenter = — (3.19) 

g 
Tracing the ray intensity is based on conservation of energy. A hypo- 
thetical amount of radiated power trapped in a ray tube produces intensity which 
is inversely proportional to the cross sectional area of the tube. Tracing the cross 
sectional area along the tube, using the ray tracing equations, provides the required 


intensity computation. 
[(D,,R1) _ S(D2,R2) 
——<—_—_ = (3.20) 
L(D2,R2) (Di ,R1) 
A singular case arises when the cross sectional area of the ray tube 
reduces to zero and drives the intensity to infinity. This is called a focal point. 
Focal points in which the intensity is very high are actually encountered in the ocean 


and the acoustic intensities measured there are very high. Ray acoustic intensity 


computations are obviously not valid under such conditions. A ray going through 
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a focal point is phase shifted by 90 deg. This phenomenon which is experimentally 
measurable, can be explained analytically although not in a trivial manner. 

Ray equations allow tracing a ray in a medium of constant SVP gra- 
dient, however, a constant SVP gradient is rarely the case in practice. When the 
SVP is not linear, it can be represented by a piecewise linear function which has 
constant gradients within depth layers. Eq. (3.14) through (3.17) are then used to 
trace a ray inside a given layer. The continuity of the linearized SVP and Snell’s 
law ensure ray angle continuity in the transition between layers. Consecutive trac- © 
ing of the ray inside the constant gradient layers can thus construct the complete 
ray path throughout the entire medium. 

c. Effect on MP 

(1) Time Delay. An example will help demonstrate the effect of 
the IH propagation on MP depth and range measurement. A set of time delays 
for an IH case is computed for a given target and set of ocean conditions using 
Eq. (3.16). These time delays are then used as independent variables of the 
homogeneous inverse function Eq. (1.2). The erroneous inverted depth and range 
are compared to the original values. 

A constant gradient of g= 0.05 sec™! in the 100 m near the 
surface,*in water depth of 513 m and a target at depth of 100 m and range of 
7 km,** produce TDOAs to a receiver at depth of 162 m of 7 = 5.48 ms and 
T = 28.92 ms. When substituted into Eq. (1.2), which is based on straight line 
propagation, a depth of 19 m and range of 5.64 km are produced. This represents 


a percent error of 20% in range relative to the true 7 km range. 


* See case C2251 in Appendix E 
** Measured relative to the receiver, see Fig. 2.1 
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(2) Lack of Closed Form Solution - The Eigenray Problem. 
Of significant importance is the fact that in all the tracing equations the initial 
transmission angle at the source (/;) is assumed to be known. This angle, which 
can alternatively be replaced by the reception angle (2), actually ‘selects’ the ray 
to be traced. 

The reverse problem of finding the terminal angles of a ray pass- 
ing through two given end points is called the eigenray problem. The eigenray 
problem does not have a closed form solution in the IH case, especially not in the 
multilayered SVP case. This is part of the difficulty which prevented compensation 
for the medium inhomogeneity in previous work witd thus became a key challenge 
of our work. 

(3) Loss of Direct Ray. Another outcome of the IH propagation 
is the possible loss of direct path between the source and receiver due to the ray 
bending. Fig. 1.3a, which was computed by a Fortran program [26] using the ray 
tracing equations, shows an example of a receiver Rx placed in an area not reached 
by direct rays from source J, which is referred to as a “shadow zone.” The first 
arrival in this case is the reflection from the bottom. The second arrival is the ray 
reflected from the surface and then from the bottom. The third arrival is the one 
reflected from the bottom first and then from the surface. The above mentioned 
difficulty to predict the exact amplitude and phase of the reflected wave is further 
enhanced for a case of double reflection such as the one shown. The polarity of 
the ACF peaks becomes less predictable. Thus, it is obvious, that the IH medium 
amplifies the problem of path identification. 

In Comparing Fig. 1.3a to Fig. 1.1, the difference in the MP 
structure is mainly that the direct path To of Fig. 1.1 does not exist in Fig. 1.3b. A 


scalar called the bounce count (BC) was devised here in an attempt to characterize 
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the MP structure. It is defined as the total number of depth extremes in the set 
of the first three resolvable arrivals (two TDOAs). 
The bounce count of the straight line propagation of Fig. 1.1 is 
2 due to the single bounces of both the surface and the bottom paths. The bounce 
count of the first three arrivals of the MP structure in Fig. 1.3b is 4 due to the 
single bounces of the first and second arrivals and the two bounces of the third 
path. The bounce count therefore distinguishes between the two MP structures. 
While ambiguous as an absolute descriptor, the bounce count 
was found to be very useful in practice primarily in identifying a change in the 
MP structure. By plotting the BC over a depth range cross section of the ocean, 
it became easy to divide the plane into regions of consistent propagation as will 
be shown later. Care should be exercised, however, in detailed analysis using the 
bounce count since radically different MP structures could produce identical bounce 
counts. 
3. Ocean Data 

The travel time of paths reflected from the bottom is obviously dependent 
on the ocean depth. Error in the assumed depth will cause errors in the inverted 
target depth and range. The CRLB for the joint estimation of ocean depth and 
target depth and range was developed by Hamilton [Ref. 19] as a composite estima- 
tion problem, with depth assumed normally distributed. The impact of a constant 
depth error is demonstrated in Chapter Four. 

4. Receiver and Delay Estimator Effects 

The effect of the limited delay resolution, and the inability to associate 

delays with paths are considered here along with an outline of a hypothetical but 


realistic form of delay estimation instrumentation. 
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a. Bandwidth and Delay Resolution 

The target broadband noise has a limited bandwidth. Further, the 
high frequency portion of the spectrum tends to attenuate in a shorter distance 
than the lower frequency portion as the sound travels through the water, due to 
absorption processes. The receiver, on the other hand is constrained by engineering 
design compromises arising primarily from the beamformer, and computational 
load associated with its sampling rates. The result is that the acoustic receiver 
bandwidth is rarely ever more then a few KHz wide. 

The limited bandwidth translates into a limited resolution between 
adjacent time lags. If one assumes a uniform signal spectral density over the entire 
frequency band of 0 to B Hz, the main lobe of the autocorrelation function of the 


original signal z(t) will be of the form 
AC F,2(7T) = ACF,,(0) - sinc(r- B) (3.21) 


and its first zero occurs at 1/B sec. Replicas of the original signal which are less 
then 1/B apart will form a nonresolvable combined ACF peak. (In the extremely 
unlikely event of two opposite sign equal amplitude replicas the ACF peak may 
even disappear altogether. ) 
b. Nonidentifiable Path 

The combination of the effects discussed earlier render the association 
of the TDOA with the corresponding path a difficult task indeed which is imprac- 
tical in many target tracking situations. The polarity of the ACF peaks is strongly 


affected by the following: 
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e Potential loss of the direct reference path. 
e Complex frequency and angle dependent reflection coefficients. 


e Phase shift through a focal point. 


An alternative approach taken here is to rely only on the time difference values, 
and to disregard the ACF polarity and its possible association with a specific path. 
The response of the time delay tracker to this nonresolvable nonasso- 
ciable situation will determine the behavior of the target tracker. For this reason, 
the specific instrumentation of the time delay estimator assumed here is briefly 
described below. 
5. Instrumentation Outline 

A discrete set of ACF values spanning the range of expected lag time, is 
computed by time or frequency domain methods, using the time sampled acoustic 
signal. A polynomial is fitted to the sampled ACF lag data and its first two 
amplitude maxima, are located (excluding t = 0). The lags times corresponding 
to these peaks are the MP time differences of arrival (or delays). The delays are 
usually further filtered by a simple exponential averager.* 

The first two lags are preferred for the inversion since they result from 
simpler paths and are expected to have higher SNR. But when any two of the 
first three ACF peaks coincide only one combined and somewhat erroneous peak is 
formed. An example where this problem is significant and very typical is a source 
placed near the surface. The direct and surface reflected paths from the source 
to a distant receiver will have very similar path length and travel time, and their 


ACF peak will not be resolvable. 


* Another technique of similar consequence uses interpolation between 
adjacent ACF samples to reconstruct the continuous delay, Friedlander [16]. 
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Other lags corresponding to more complex and longer paths, do exist and 
can be used in such cases. Having longer lag time, these delays will not vanish, due 
to loss of resolution, incidentally with the first two. The third or higher lags can 
thus provide an alternative measurement, at least while one of the first two lags is 
lost. 

In the instrumentation model simulated here a third ACF lag is tracked all 
along, and is substituted as the second lag only when any of the first three arrivals 
coincide. The two nonresolvable lags are reported as one lag under this situation. 
Coincidence is defined with some safety margin that ensures substitution of the 
third lag prior to the actual loss of resolution. 

In Fig. 3.3 plot of the first three TDOA as a function of history time k 
is shown. At time Ty, the ACF peak of 7% is not resolvable from the main lobe 
at rt = 0, and is consequently replaced by 7) which in turn is replaced by 73. The 
third TDOA 73 was tracked but not output before Ty2. If 7 and 7 merge at Ty, 
as shown in Fig. 3.4 their combined peak will be reported as 7 and 73 will replace 
7. 

In the straight line homogeneous propagation, identification of the mea- 
sured delays as 7 and 72 is assumed such that they can be used in Eq. (1.2). In the 
realistic IH case the delays are identified instead as the first (shorter) and second 
(longer) lags. A vector of two time delays, the first of which is always the smaller, 
is output by the estimator and these modified TDOAs are denoted as 7; and 7%. 

6. Resulting Ambiguity 

The proposed method of time delay extraction immediately raises the ques- 

tion of ambiguity. There can obviously be more then one position with the same 


two time delays as demonstrated by the simple following case. 
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For an observer at depth 162 m and water depth 503 m the targets: 


Position Target 1 Target 2 
Depth m -100 -242 
Range m 6000 4782 


produce 7o,7,,72 direct surface and bottom travel time, respectively of 


Travel time [sec] Target 1 Target 2 

To 4.00055 3.192079 

j i 4.00997 3.210297 
T; 4.01877 3.201499 


yielding travel time differences of 


TDOA [msec] Target 1 Target 2 
™| =1T; —To 9.41 18.2 
Q= T> — To 18.2 9.41 


After sorting these become 


Modified TDOA [msec] Target 1 Target 2 
ti = Min{n, 7} 9.41 9.91 
to = Maz{n, 72} 18.2 18.2 


Yielding the same 772 set for two vastly different depth range positions. 

This ambiguity is a direct result of the inability to associate paths and 
delays. Means to reduce this ambiguity to a manageable proportion are covered in 
Section C. A more effective use of the high bounce count paths designed to further 


reduce the ambiguity is discussed in Chapter Six. 
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7. Graphical Representation of the Medium and Receiver 
Effects 


a. Introduction 

The 2-D direct function transformation of source depth and range 
(DR) position, to MP time delays, (7,, 72), can be viewed as surfaces of delays over 
the 2-D depth-range position plane. These surfaces are the basis of the MP inver- 
sion procedure and understanding their nature is important in the development of 
the inversion process. 

The graphical interpretation of the direct function will be presented 
in this section including the following effects which the inedium and receiver have 


upon this function: 
e Lack of path identification. 
e Limited resolution. 


e IH propagation. 


Four types of surfaces are presented starting with a theoretical case of 
homogeneous propagation, and resolvable and identifiable multipaths. The above 
listed effects are then added in three accumulating steps ending with the realistic 
case of IH propagation, limited resolution and non identifiable paths. 

b. Case Definition 

The number of parameters determining the MP structure is very large 
and attempting to investigate all of their various combinations is an impossible 
and not a valuable effort. A number of representative cases were selected and 
the parameters involved were coded into cases. A case name is composed of one 
letter and four digits, e.g., C1123, which has the following meaning. The leading 


character defines the type of assumptions and instrumentation used. The first 
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_ digit defines the SVP and the second digit defines the ocean and observer depth. 
The last two digits define a grid of depth range points for which the TDOAs are 
computed. A detailed description and settings of the various parameters is given 
in Appendix E. . 

To demonstrate the nature of the direct function surfaces three cases 
were used. The first case for an ocean depth around 500 m and receiver depth of 162 
m (U1111) represents homogeneous propagation and infintely resolvable TDOAs 
which are associable with the paths for bottom depth of 503 m and receiver depth 
of 162 m. The second case (51111) represents the same case but without the 
associability of TDOA with the multipath. The third case (C1111) represents 
the same case with both the TDOA-path associability and the perfect resolution 
assumptions removed. Three dimensional (3-D) surface and contour plots of TDOA 
1,72 or t},t2 over the DR plane are plotted for each case. Regions of consistent 
MP structure are indicated by contours of bounce count over the DR plane. A 
short description points to the issues of interest in the plots. 

c. TDOA Surfaces 

Homogeneous, perfectly resolvable and identifiable U1111. 
TDOA surfaces and contour lines for this ideal case are shown in Fig. 3.5. 
Both surfaces are monotonic and 7, can be either smaller or larger then 7. Both 
delays decrease with range but opposite changes in 7 T) result when depth is chang- 
ing along a fixed range. Note that 7, is small and 7 is large close to the surface, 
and the reverse is true near the bottom. No bounce count plots are presented since 
the count is 2 all over the plane. Note the vanishing 7 and 7 near the surface 
and the bottom respectively which causes the merging of the corresponding ACF 
peaks with the main ACF lobe at r = 0. These lags are substituted by those from 


multiple reflection paths in the limited resolution case. Also note that the TDOA 
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Fig. 3.5. Ideal TDOA surface (U1111) 
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_ vanishes with range, this tends to make the effect of delay estimation noise more 
significant at longer ranges. Since in addition the TDOA estimation noise increases 
with range, the MP method becomes inherently limited to short ranges. 

Homogeneous sorted perfect resolution $1111 modified TDOA 
™,7 surfaces and contours for this ideal case where association of the path and 
delays is not assumed are shown in Fig. 3.6. When the TDOAs are sorted to 
produce t, and ¢2 as described earlier the result is tf; < t2 by definition. The 
monotonocity of the surfaces is lost along the DR loci where t, = t2. The trend of 
reduced delays 7,72 at long ranges has not changed. The depth still determines 
the relation between 7 and 7) but not in the simple unambiguous manner as in 
the unsorted case. 

It is only 7, now that vanishes near the surface, when it becomes 
smaller than the delay estimation resolution it will be replaced by what is here 
plotted as 7. The 7, will then be replaced by the arrival lag of the next multipath, 
the one bouncing from both the surface and the bottom. Since this will only happen 
for the realistic receiver, the bounce count here (perfect resolution) is still 2 all over 
the DR plane and is therefore not plotted. Contours of TDOA in the first 100 m 
near the surface are shown in Fig. 3.7 as a reference for the later comparison to an 
IH case. 

Homogeneous realistic finite resolution C1111. Fig. 3.8 shows 
the resulting surfaces when the TDOA are not associable with the multipath 
(sorted) and a finite delay resolution (0.5 ms) is assumed. 

The main difference between this and $1111 case is the substitution 
of a high order reflection time delay when time delays from simpler reflections are 
not resolvable. This comes about near the surface where 7, is small and at depth 


of around -176 m relative to the receiver where 7 and 72 are similar causing the 
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_ replacement of 72 by the lag corresponding to the higher bounce count path (#3). 
The bounce count plot is helpful in recognizing this effect which of course is also 


obvious on the TDOA surface and contours as well. 


C. IMPROVED FILTER 
1. Concept Overview 

The subject of this section is inversion of realistic IH multipath time differ- 
ence of arrivals (TDOA) to source depth and range coordinates. Since there is no 
closed form solution for either the direct or the inverse relation between time delays 
and position, a two-step numerical solution is developed. As before the function 
that computes depth and range from time delays is referred to simply as the direct 
function, and the procedure that computes depth and range from time delays is 
referred to as the inverse function. The latter may at times be multivalued and so 
may not be a one-to-one function in the mathematical sense. 

First, the direct function values of the TDOA are computed over a dis- 
crete grid of target depth and range positions. This computation is done off-line. 
Measured time delays are then inverted to obtain position (depth and range). This 
computation must be done on-line. 

In the off line computation an eigenray model is used to find the travel 
times of the multiple paths between the receiver and a given depth and range 
point. The travel time produced is further processed to simulate the realistic delay 
estimation. Generation and processing of the TDOA is repeated for every point in 
the entire DR grid. The processed TDOA data thus produced, is stored as a table 


of time delay as a function of the depth and range. 


1Z1 


In the on line step, a set of two measured time delays is inverted to produce 
depth (and range) using 2-D linear interpolation. Depth (and range) values of 
three points in the direct function table, which have time delay values closest to 
the measured set, are used as a basis for the interpolation. To locate these points 
a unique search algorithm is developed. The algorithm takes advantage of a priori 
information available in the form of the predicted target position. The prediction, 
in turn, is derived from past depth and range measurements processed by the target 
tracker. 

The use of the eigenray model and the delay estimation simulation to 
generate the direct function is described in Section C.2. The inversion process is 
described in Section C.4. 

2. Direct Function 

The precomputed tabulated direct function of the modified TDOA over 
the DR grid is computed by means of an eigenray model which is repetitively 
exercised over the DR grid points. This process is described next. 

a. The Eigenray Model 

An eigenray model is a program that finds the multiple rays which path 
between two given end points, providing their angles, path length and travel time. 
Computing the TDOA of the MP structure requires very accurate computation 
of the path travel time. This results from the TDOA being a difference of a few 
milliseconds between travel times which are of the order of a few seconds. 

In the search for an eigenray model only compact programs requiring 
small computation time on readily available microcomputers were considered. This 
was done in order to maintain compatibility both with the computational resources 
available to this research and those available in onboard systems. Since no closed 


form solution for the eigenray problem is available (recall discussion in Section 
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III.B.2.c) the eigenray models are iterative in nature. The conflicting requirements 
for high accuracy, limited word length and short run times are met by only a few 
existing models. 

Appendix F discusses the SMART ( SMall Acoustic Ray Trace, Ref. 
22) selected as the eigenray model for this thesis*. The program computes the 
dependence of range on initial angle and uses this relation in its iterative search 
for the eigenray. 

b. Generation of Direct Function Table 

The SMART model produces the sound travel times for an eigenray 

set pertaining to one receiver and one source. The generation of the complete direct 


function table requires the following five additional steps. 


1. Repeatedly exercising the model for all the source points on the discrete grid 
of target depth and range (DR) positions. 


2. Selecting the shortest travel time as the first arrival and subtracting it from 
the travel time of the other paths. The resulting differences correspond to the 


time lags of the ACF. 


3. Sorting the paths in ascending travel time order. This corresponds to the 
estimator logic which assumes that it is not possible to directly associate the 


TDOAs with the paths. 


4. Applying the resolution limitation by eliminating delays which are less then 
the resolvable difference away from their preceding delay. 


5. Selecting the shortest two delays and storing them as the tabulated direct 
function. 

The last step is in need of further elaboration. The data structure used 

for the purpose of storing the direct function table is a 3-D array Q of dimension 

2x Na x N; indexed by @,2, 7. The first array index @, selects the first and second 


* The proprietary model was made available to this research courtesy of 
Mission Sciences Corp., Commack, N.Y. 
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delays t,;t2(@ = 1,2). The second and third indices 2 and j select a point on the 
depth and range grid respectively. The variables Ng and N, are the number of 
depth and range steps in the grid. Each point [z 7] in the grid thus conceptually 
corresponds to a four element vector called a quad (q) which has the following four 


values: 


qij = D,|2], Ry[z], ti (2, 7], to(z, 7] (3.22) 


where D,[z] is the depth implied by the depth grid index 2. R,[j] is the range 
implied by the range grid index 7. t,, t2 are the two modified TDOAs given by ¢; = 
T, —Tp, and t2 = T,—T7p, and stored in positions Q[1, z, 7] and Q[2, z, 7] --spectively, 
where To, T;, Z2 are the ordered travel times of the first three resolvable paths such 
that Tp < ZT, < 72. Namely 


tie toe (3.23) 


The vectors D, and R, are the depth and range grid scale vectors 


D, = Dy + AD-(0,1,2,.. ea (3.24a) 


R, = Rg + AR-(0,1,2,...Ne—1 (3.246) 


with Doo, Rgo the initial grid points at shallow depth and short range, and AD, AR 
are the grid step sizes. 
c. Interpolation 
Linear interpolation of the TDOA values for points of the DR grid is 
used to provide a continues TDOA function of the DR position. The values of the 


t,(t2) depth and range of the three quads 


Gi, 115 Vig Gite J8 
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which form the closest triangle that surrounds an input point D;, Rin; on the DR 


plane, are used to linearly interpolate the ¢,(t2) values at that point. Two linear 


planes of the form 


T=AD+BeR+C, 


are formed in the t; DR and t2DR subspaces as follows. 


Teltinn| = Ae- Di + Bei + Ce 


T¢[t2J2|] = AeD2 + Bee + Ce 


T¢[t3)3] = AeD3 + BeRz + Ce. 


If one defines the vector Ty as 


Te = Te|t1 ji], Te[22J2], Te[23.3 | 


where = 1,2, then Eq. (3.26) can be written in matrix form as 


D, Rf 
Te=|Do Fo 
D3; Rs 





Solving for the parameters one obtains 


Ag D, 
Be | = | D2 
C'e D3; 





| 


Ay 
Bigaln. 


Ce 


(3.25) 


(3.26a) 
(3.266) 


(3.26c) 


(3.27) 


(3.28) 


(3.29) 


The plane parameters the t,t2 at any specific point D;,R;, can be interpolated 


with the formula 
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4 =(Ay By Cy) es (3.30a) 


and 


te =[A2 Bo C2)-[Din Rin ie (3.305) 


as shown graphically in Fig. 3.9. 
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Fig. 3.9. Linear interpolation. 


The ordering of the tauuiated discrete direct function according to 
depth and range enables finding the triangle of closest quads by rounding the 
input D;, Rin to the closest DR grid values. 
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d. Implementation 
SMART is a proprietary product which was obtained for this research 
in a form executable on an IBM personal computer. The travel time data generated 
by this program was transferred to a large mainframe computer (IBM). Here, the 
rest of the processing in steps 2 through 5 and the target tracking simulations 
were performed. The mainframe was selected primarily because of its good APL 
interpreter and the excellent engineering graphics support provided by GRAFSTAT 
[Ref. 27]. 
3. TDOA of IH Surfaces 

Two examples of the resulting interpolated discrete direct function surfaces 
are presented here. The form is the same as that used for the homogeneous cases 
presented in Section B. The cases are coded using the casename code described in 
Appendix E. 

An inhomogeneous finite resolution nonassociable case (A2251) is pre- 
sented first. The TDOA surfaces of SVP with a positive surface gradient of 0.05 
sec”? and very high delay resolution (0.05 ms) are shown in Fig. 3.10. Note that 
the DR grid only covers the first 100 m of the 500 m water column. The contour 
plot of the ¢; surface for the homogeneous case (case $2251) was shown in Fig. 3.7. 

While the basic trends of the ¢,¢2 surfaces are similar for homogeneous and 
IH cases, the specific TDOA values differ. This can be seen most distinctively on the 
contour plot for t; which here happens to be the delay between the surface bounce 
and the direct path, and is in general longer than the t; produced by straight line 
propagation. The difference results from the fact that the average speed of sound 
is smaller along the IH surface bounce path than along the homogeneous surface 
bounce. The bottom bounce lag ?2 is less affected since the most of its path is in 


the region below the observer where the speed of sound is constant and identical 
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Fig. 3.10. Inhomogeneous TDOA surface (A2251). 
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in both cases. The TDOA surfaces for the same medium conditions, but with a 
finite delay resolution of 0.5 ms (Case 2251), are shown in Fig. 3.11. The inability 
to resolve the surface bounce is clearly seen along the first depth grid line of both 
TDOA plots. 
4, 2-D Function Inversion 
a. Overview 

The inversion algorithm is an on-line transformation of a 2-D TDOA 
input measurement vector into a 2-D depth and range position vector DR. The 
time delays need not, and in most cases do not, coincide with points on the direct 
function generating grid. 

The task can be stated graphically when the 2-D transformations are 
viewed as mappings between the position and delay planes. The direct function 
maps constant depth and range lines (DR grid) into curved contours on the t;t2 
plane; the inverse function maps a t,¢2 grid into curved contours on the depth 
range plane. The inversion task is thus defined as finding the intersection of the 
inverse mapping ¢; tz contours and reading its DR coordinates. 

Typical inverse mappings of constant t,t2 on the DR plane were shown 
in Fig. 3.12 graphically demonstrating the inversion process. Two contours corre- 
sponding to t; value of 7.5 ms and f2 values of 30 and 50 ms are shown on Fig. 
3.12. One intersection of the 7.5 and 30 ms contours is shown (a). Four intersection 
points of the contours corresponding to t,t2 values of 7.5 ms and 50 ms are shown 
(b through e). The multiple intersections indicate the inherent ambiguity of the 
realistic case which was discussed in Sec B. 

Typical direct mappings of constant depth and range (DR grid) into 
the ¢;t2 plane are shown in Fig. 3.13a. These contours are computed using 


ideal straight line propagation. From Eq. (2.106) it is clear that the depth is 
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Fig. 3.12. Graphical inversion. 
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dependent on the ratio 7/72 while the range is dependent on the weighted sum 
(D, — Do)m + Dot. This dependence explains the straight line nature of the 
inverse direct mapping. The dots along the contours of Fig. 3.13 represent quads. 
The depth and range are defined by the contours and the TDOAs are the Cartesian 
coordinates of the dots. 

A typical direct mapping of quads for an IH realistic case is shown in 
Fig. 3.13b. Note that the order of the points on the plane is lost and that the 
points are restricted to the half plane defined by the t2 > t,. This results from the 
definition of the modified TDOA (Recall Section B). 

A cornerstone of the inversion is the fact that a quad q;; 1s in itself 


invertible in the sense that it provides the transformation 


ty } (2 +) 
= 3. Hl 
aay 7S a _ 
as well as the inverse relation 


(au) ~ ( ea) | (3.32) 


The quads can therefore be viewed in reverse, namely as defining depth and range 
values over the t,t2 delay plane. Linear interpolation can be called upon again 
to compute depth and range values at ¢;t2 points which do not coincide with any 
quad. 

Although the ‘nversion procedure seems theoretically quite simple, its 
practical implementation is complicated by the following fact: the quads in the Q 
array are not ordered on the t;t2 plane as they are on the DR plane grid. This 
is especially true for the realistic [IH case as — in Fig. 3.13b. The search for 
the three interpolation quads is therefore much more elaborate than the simple 


rounding procedure used for the direct interpolation. The retrieval of the quads 
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Fig. 3.13. Depth range mappings. 
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here involves a search which is drastically effected by the ambiguity of the surfaces. 
The search is described in the next section. 

Once an appropriate quad triangle is located, depth and range are 
interpolated. If multiple solutions are found one has to be chosen and reported 
as an output. The complete procedure includes the search, the interpolation, and 
the ambiguity resolution. The inversion is defined here as the inversion procedure 
(which is performed by the prefilter). The inversion procedure is described in 
Section C.4. 

b. TDOA Search 

Every three direct function quads, each relating values of depth D (or 
range FR) to delays t,t, can specify a plane in the Dt,t2 (or Rt,t2) space which 
can be used for the inverse interpolation. In order to improve the interpolation 
accuracy, only quads forming small triangles which surround the input t,t, point 
are considered. The total number of combinations of three quads out of the total 
(ND x NR) quads is (NPN) This makes exhaustive search for the smallest 
surrounding triangle an impossibility. 

Since different triangles produce different interpolation results, succes- 
sive prefilter iterations for the same t,t2 point, may produce different depth and 
range values if the selected triangles are not the same. This is an undesirable situ- 
ation. For example, a static target could produce changing position measurements, 
causing the target tracker to respond as if the target were maneuvering. 

The impact of a nonunique triangle selection is demonstrated in Fig. 
3.14. Part c and e of the figure show the same four quads (q) 429394) grouped in 
different triangles together with the resulting interpolated range. The quads are 


placed on a cartesian grid for clarity purposes only and the values marked near the 
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Fig. 3.14. Non-unique triangle selection. 


quads (8 and 12) are the range values of the quads. Three interpolation triangles 
can be selected for range interpolation for a target moving from point e to point c. 


These triangles are the set (case A) 


big) (q3, q2, 94) 
rte (q1 143; q2 ) 


or the single triangle (case B) 


Iri3 (q1 1 93; q4) P 


Selecting T'rz1 for target position e to d and tri2 for position d to c (case A) will 


provide a continuous range interpolation. Selecting trz3 for the complete trajectory 
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a to c will also yield continuous interpolation (a constant). Selecting T'rz3 for posi- 
tions e to d followed by trz2 for d to c will result in the discontinuous interpolated 
range shown in d. This discontinuity in range is troublesome for the maneuvering 
target tracker. 

The uniqueness problem has two aspects. One part of it results from 
the unordered nature of the points on the t,t, plane. This part can be overcome 
by constructing unique triangles either for the whole grid at one time or locally, 
when and where needed using some algorithm. An example of this problem is the 
fact that every triangle that surrounds the input point can itself be surrounded by 
another yet larger triangle. A simple rule of choosing the smallest possible triangles 
(for example by area) could solve this problem. 

The other aspect of the problem, is the inherent ambiguity resulting 
from the non-monotonicity of the direct function. This ambiguity, discussed earlier, 
is demonstrated here again with an emphasis on its relation to the triangle selection. 
Six quads marked gq, through gg are shown mapped on the ¢;t2 plane in Fig. 3.15a 
with the corresponding depth values as the vertical axis. The folding surfaces 
formed by the quads helps one understand the name Manifold by which they are 
called. As expected, the non-monotonicity of the direct function results in two 
possible interpolation triangles trzl: (q1,q2,q3) and trz2: (q4,q¢5,q¢6) for an input 
point-p. The problem arises because both triangles are “valid.” Any attempt to 
resolve the dilemma by constructing a “unique” triangular grid like the sample in 
Fig. 3.15b may create wrong interpolating triangles by mixing delay values from 
two separate areas of the non-monotonic function (e,g. 42, 44, 43). 

The two triangles represent two actuate areas on the DR plane for 
which the particular MP structure generates similar delay values t,t2. This is 


caused by the following combined effects: the IH propagation; the limited resolution 
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receiver; and primarily the lack of path delay association and the corresponding 
assumed instrumentation logic (described in Section B). The triangles selection 
mechanism that is used, has to retain the option of selecting both triangles tr7z] 
and trz2 since they are both valid for interpolation. 

Attention is directed at this point to the fact that while trz1 and 
tri2 are partially overlapping in the TDOA plane, they are disjoint in the DR 
plane. This difference is significant and is used later to assist in the selection of 
the triangles. 

An algorithmic solution to the general problem of uniquely selecting 
a triangular grid over a set of unordered points in a plane does exist. Surprisingly 
enough, it was developed only recently (1975) by Shamos [Ref. 28]. The solution is 
elaborate and involves a geometric construction but it does guarantee unique trian- 
gles. However it does not solve the above mentioned problems of eliminating valid 
and creating nonvalid triangles. The Shamos solution is therefore not applicable 
for the problem of non-monotonic function inversion and a different new approach 
had to be devised here. This approach is described next. 

c. DR Search 

The difficulties encountered in selecting triangles resulted from the 
unordered scattering of the quads in the TDOA plane (Fig. 3.13b). On the DR 
plane the quads are ordered on the Cartesian grid. Constructing a unique disjoint 
set of triangles which exhaustively covers the depth range plane 's trivial and can 
be done by splitting each rectangle of the DR grid into two triangles using parallel 
diagonal lines as shown on Fig. 3.16. The set of triangles thus formed is defined 


as the DR triangles space. A key idea in the search is to consider only triangles 
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Fig. 3.16. Unique valid DR triangle set. 


in the delay plane that arise from the mapping of triangles from the DR triangles 


space. This delay triangles space has the following distinct advantages: 


si 


2. 


The set 1s unique. 


The set is completely valid. That is, the triangles are locally the smallest and 
will not be constructed by erroneously grouping quads from disjoint areas of 
the direct function. 


. The set is sufficient. That is, the option of choosing all valid TDOA triangles 


is Maintained. 


. The search space size is significantly reduced. The number of possible triangles 


is reduced from (Ne Ne) to 2-Nq-N,. This is a reduction from order O(n!) 
to order O(n). 


. Construction of the set is trivial and does not require computation. 


. The set facilitates an efficient structured search. The search for the surround- 


ing TDOA triangle can utilize the structure of the underlying DR grid in an 
effective binary type search. 
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Viewed in the context of the graphical statement of the inversion prob- 
lem (Fig. 3.12) the search now is for the DR triangle that surrounds the intersection 
point of the t,,¢2 contours of the inverse function on the DR plane. The existence 
of a solution is conditioned by the monotonicity of the direct function over the 
region of the triangle which dictates that the DR grid be made smaller than the 
minimum monotonic region. 

Since the DR grid of quads is Cartesian, searching for a minimum rect- 
angle that contains the desired triangle is more convenient. After the surrounding 
rectangle is found the specific triangle(s) can then be located by examining the two 
triangles that make up the rectangle. 

d. Localized Search 

Finding the particular triangle in the total set of 2-Nq-N, grid points 
that surround the ¢,;t2 input measurement point requires further searching. The 
test to see if the point is within a triangle,* involves some computation and the 
search space, though reduced can be still very large. Hence the possibility of ex- 
haustively considering all the triangles is not desirable. A binary-type search made 
possible by the DR ordered grid structure can reduce the average computation 
load. In addition it can focus the attention on the area where the solution is most 
likely to be found. 

A place at which to start the search (a root) is provided by the target 
tracker predicted position. If a target state vector has already been established, 
then the next target measurement is most likely to be found in the vicinity of the 
predicted DR position. Starting the search with the rectangle that surrounds the 


predicted position is only the natural choice. 


* The test, referred to as the inversion test is described later. 
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The Q array indexes required to retrieve the quads associated with 
this rectangle are produced, as they were for the direct function interpolation, 
by rounding to the closest DR grid values. If the rectangle that surrounds the 
predicted point also surrounds the t,t2 point the search has produced the required 
triangle and it terminates. Upon failure, however, the search has to continue over 
a growing and more distant set of DR rectangles until a surrounding one is located. 
This local search ensures that the solution closest to the target will be found first. 
This anifcanitl} reduces the ambiguity inherent in the MP inversion. 

Two different search strategies were designed and implemented for this 
purpose. The first, titled MEA (means ends analysis) is a directed search strategy 
common in artificial intelligence applications. The second method referred to as 
local search algorithm (LSA) is a breadth first type search performed in tiers of 
growing distance from the predicted point. 

The MEA method performed very effectively for monotonic direct 
functions but failed in non-monotonic cases. The LSA algorithm performed very 
well for all cases and the computational load, usually associated with breadth first 
search, was reduced by an innovative dual phase evaluation algorithm. A brief 
description of the local search algorithm is presented next and a more detailed ex- 
planation of its components is presented in Appendix G. The system block diagram 
is shown in Fig. 3.17. Note the feedback of the predicted position as a root for the 
prefilter’s local search. 

e. The Local Search Algorithm (LSA) 

Given the TDOA measurement point the goal of the LSA is to produce 
an interpolation triangle that surrounds the intersection of the inverse function t,t 
contour lines. The search commences in a DR rectangle where the solution is most 


likely to be found, namely at the predicted target position. The search progresses 
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Fig. 3.17. System block diagram. 


oy increasing the rectangular search area until a solution is found or until the 
dimensions of the area have reached a specified maximum limit. Upon termination 
(success or failure) the algorithm reports the number of solutions found, the size 
of the area actually searched and the triangles found (if any). 

The localized search is a part of the overall inversion algorithm. The 
inversion initiates and uses the outputs of the search to produce the desired mea- 
surement depth and range. The interface between the inversion and the search is 
given in Table 3.1. 

A numerical analysis procedure, originally developed for finding zeros 
of complex functions (Hamming [29]) was modified to apply to the search. One test 
and two recursive procedures are used in the LSA, they are the inversion test the 
outward and the inward searches. The inversion test examines a given rectangular 


on the DR grid to find out if a solution may exists in it. The outward search 
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TABLE 3.1 


LSA INTERFACE VARIABLES 


LSA input: 
tite Delay measurements (input point). 

DR prd. Target predicted position, depth range. 
Max size Maximum search area size ,depth range. 
LSA output: 

Num The number of triangles found. 
Area size The size of area searched before 


the solution was found. 
Solution The quad triangles. 


examines the current (input) rectangle using the inversion test. As long as the test 
fails (no solution in the current area) and as long as the maximum search area has 
not been reached, the search area grows outward. If the maximum area is reached 
and a solution was not found a failure is reported. 

The inward search attempts to locate a solution inside a given rect- 
angular area. It starts by performing the conditional inversion test on the entire 
area. If the test succeeds the area is divided to four quadrants and a solution is 
searched for in each quadrant by means of a recursive call to the inward search. 
The recursion terminates either when a 1 x 1 rectangle is reached, or when the 
test indicates that a quadrant does not contain a solution. If a 1 x 1 rectangle is 
formed, an additional surrounding test is performed on the two DR grid triangles 
that comprise that rectangle to determine which if any of them contains a solution. 
All the triangles thus found are reported as outputs to the inversion procedure for 
further processing. 

The search as a whole is thus completed. The surrounding triangles 


closest to the predicted target position are reported to the inversion procedure 
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for the further processing which is discussed next. The local search algorithm 
is demonstrated in Fig. 3.18. Fig. 3.18(a) shows two iterations of the outward 
search initiated at point p. Fig. 3.18(b) shows three iterations of the inward which 
produce two solutions corresponding to the two intersections of the t,t, contour 
lines. A more detailed description of the LSA is presented in Appendix G including 
a detailed pseudo code description. 
f. The Inversion 

The inversion procedure performs the prefilter operation of transform- 
ing the input measurements TDOAs to depth and range measurements. It uses the 
direct function table and the local search algorithm described earlier. The inversion 


performs the following five tasks. 


e Initialization of the search 


Interpolation 


Ambiguity resolution 
e Failure handling 


e Performance monitoring. 


Each of these is described in more detail below. 

(1) Initiation of the Search Process. The process is initiated by 
setting up the first 1 x 1 search rectangle and starting the outward search. The 
initial rectangle is set by quantizing the target predicted position to the closest DR 
grid position. 

(2) Interpolation. Once the appropriate three quads are located. 
the depth and range need to be interpolated. Lise and bilinear interpolations 


were used. 
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Fig. 3.18. Outward and inward search. 
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Linear interpolation is done in much the same way as the time 
delay interpolation is done. The interpolation equations are similar to those de- 
scribed earlier for the direct function interpolation. (See Section C.2.c.) 


A linear plane in the depth (range) t,,t2 space of the form 


D = Aat; + Bat2 + Ca (3.33) 


is defined by the three quads q721, 1; g?2J23 q%3J3 Where q;; is defined in Eq. (3.22) 
and the subscripts index the three points of the selected triangle on the original 
depth range grid. The equations for the plane coefficients in the depth t,t) space 


are 


Dita ji] = Aa - ty [2191] + Ba - te[tigi] + Ca 
D{t2j2] = Aa: ti[t2j2] + Ba: to[tejo] + Ca (3.34) 


D{t333] = Aa: tilt393] + Ba - te[¢373] + Ca 
which takes the matrix form 


Ag ty[tigi] teltigi] 1 Dlti91] 
Be| = | t:féeg.] tala. 0 - eel (3.35) 
C4 ti[?373] te[23j73] 1 D{t3 33] 


from which the parameter vector (AqBygC a)? can be solved. The depth (range) at 


the specific t)_t2_ measurement point can be interpolated as 
D=(AgBaCa)-(t1 t2 1). (3.36a) 
And similarly for range 


Rox (Ap De Co): (Gare. “er (3.366) 


with A,B,C, solved by an equation like Eq. (3.35) with R[z,7] substituted for 
Diego |. 
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Bilinear interpolation was tested as well. The bilinear interpo- 


lating surface [Ref. 30] applied to depth interpolation takes the form 


D=Aq:ti + Ba: te +Ca-t, -t2 + Dy. (oe3i1)) 


The coefficients Ag, By,Ca, Da of the bilinear interpolation can be determined by 


solving the matrix equation 


ty[tigi]’ toftagi] ta(21gi]-teleizi] 1 Ag D213, ] 

ti[t2j2] te[t2j2]  tilt2g2]-telteze] 1 Ba} _ | Dit232] (3.38) 
ti[2373] telzsj3] ti{tsj3]-te[e373] 1 Ca| | Dlts33] 
ty[taja] tof[taja] tilts]: toltaza] 1 Da Dita ja] 


While linear interpolation uses a triangle formed by three quads, the bilinear inter- 
polation requires a square of four quads. The closest square, however, is naturally 
produced by the local search algorithm. 

(3) Ambiguity Resolution of Multiple Solutions. When multi- 
ple solutions exist in the search rectangle a choice is made among them by interpo- 
lating all of the candidate solutions and selecting the one closest to the predicted 
position (point $2 in Fig. 3.18). For a single linear measurement and a normally 
distributed measurement noise and prediction error, this solution will be the more 
likely one. Realistically the measurement is nonlinear, and the two measurements 
(depth and range) do not have the same variance. An alternative more rigorous 
choice would imply extensive computations that are expected to produce practically 
the same results. A potentially more rigorous approach to the overall ambiguity 
problem is discussed in Chapter Six. 

(4) Failure Handling. A search can terminate without a solution 


for one of the following two reasons: 


1. The measurement noise resulted in placing the t; t2 input point outside of the 
maximum search area or even outside the complete DR grid. 
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2. The search area reached the maximum size without covering the actual mea- 
sured source position. 


When an inversion triangle is not found, the predicted target 
position is returned by the inversion procedure as an output. This can be thought 
of as a lack of innovation in the current measurement, and effectively reduces the 
sampling rate of the system. The tracker performance degrades as a result. When 
the search failure occurs infrequently (low failure rate) only the response to fast 
target maneuvers is hampered. Frequent failures (high failure rate) can lower the 
effective sampling rate of the observer below the Nyquist rate required to sample the 
target’s dynamics. This drastically reduces the quality of the tracking. Frequent 
failures are unfortunately fed back in the form of search areas centered around the 
wrong target position. This in turn increases the chances of failure and eventually 
leads to a complete loss of tracking. A large maximum search area limit will reduce 
this effect, but requires increased on-line computation. The maximum search area 
limit is thus a key system parameter which needs to be optimized for the prevailing 
noise conditions and the available computational resources. 

In experimental simulation of the MP tracking system the overall 
effect of the inversion failure was surprisingly gradual and meaningful tracking of 
maneuvering targets was demonstrated with up to a 50% failure rate. The robust 
performance can probably be attributed to a more-or-less symmetric delay error 
distribution. In such cases the elimination of the extreme measurement errczs 
at either ends of the distribution caused by the failure does not affect the mean 
depth and range values. Therefore, as long as the tracking errors are not large and 
the search area is centered properly around the predicted position, the impact of 


eliminating extreme delay measurements is mild. 
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(5) Prefilter Performance Monitoring. The inversion proce- 
dures also monitors the performance of the search by averaging the number of 
solutions and the search area size. The averaging is done by means of an auto- 
regressive lowpass filter. The average solution count should ideally be one, a lower 
value indicates search failures, while a higher count indicates ambiguity. The search 
area average helps optimize the maximum search area and is a measure of the load 
associated with the search. It also provides a good indication of the effects that 
delay measurement noise has on the search. The solution count and the area size 
proved to be very useful in the overall system performance evaluation. 

g. 2-D Function Inversion Summary 

A 2-D function inversion algorithm was devised, capable of translating 
realistically estimated MP TDOA generated by an IH medium to source depth 
and range position. An eigenray acoustic model identifies and precomputes the 
the travel time of rays to the receiver from hypothetical targets placed at the 
coordinates on a depth range grid. The travel time is transformed to TDOA using 
an assumed, ACF based, realistic receiver and delay estimation instrumentation 
and stored in a tabulated direct function. The inversion is performed by inverse 
interpolation over three or four relevant points in the grid located by means of a 
special local search algorithm (LSA). A total of three prefilters were programmed. 
They were numbered with prefilter numbers (PFN) 4 through 6 (see Chapter Two, 
Section E). Prefilters PF'N-4 used linear interpolation and the MEA search; PFN-5 
used linear interpolation and the LSA. PFN-6 used bilinear interpolation and the 
LSA. 

The procedure was carefully evaluated and the detailed evaluation 


scheme and results are presented in the next chapter. In general the performance 
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was good indicating that the procedure is suitable for incorporation into a multi- 


path tracking scheme. 
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IV. IMPROVED PREFILTER PERFORMANCE 


An intensive evaluation of the prefilter was performed and the results are 
presented in this chapter. Section A discusses the motivation for the evaluation 
and some of the expected results. Section B presents the evaluation scheme. Section 


C contains a detailed presentation and analysis of the results. 


a 


A. MOTIVATION 

The prefilter forms a central part of the overall MP target tracking system. 
Examining its performance is a necessary prerequisite to analysis of the overall 
system. The simulation has to provide all the realistic conditions under which the 
prefilter is expected to operate as a part of the MP tracker. In addition, the sim- 
ulation conditions have to be separately controlled in order to isolate their impact 
on the prefilter. The effects of the following issues on the prefilter performance are 


investigated: 
¢ Medium effects. 
o IH propagation. 
o Propagation loss. 
o Bottom depth errors. 
e Delay estimation. 
o Delay noise bias. 
o Noise range dependence. 
e Prefilter design 
o Grid size. 


o Interpolation errors. 


lol 


The overall performance of the prefilter is expected to vary over the depth and 
range grid. This is due to the range dependence of the noise and the nonlinearity 
of both the inversion and the time delay estimation algorithms. To separate the 
effects the evaluation is performed over the entire 2-D grid in two steps. In the first 
step a range independent TDOA noise is assumed; in the second step the noise is 


range dependent. 


The above nonlinearities give rise to two range bias errors. The first is a range 
underestimation bias which results from the overestimated TDOAs produced by 
the nonnegative effect (see Appendix D.) The other is caused by t!-e delay to 
position transformation, which is inherently nonlinear. This second bias tends to 
overestimate the range, and is discussed next. The effect of both biases on depth 
is dependent on the particular MP structure. 

An example of the effect of the nonlinear TDOA to range inverse transforma- 
tion, is shown in Fig. 4.1. A section of the range dependence on ft, for a constant 
value of t2 1s shown in Fig. 4.la. A zero mean normal input noise distribution 
(Fig. 4.1c) is transformed to an asymmetric and non zero mean range distribution 
(Fig. 4.1b). The typically negative value of oF leads to an overestimated range. 
This can be seen in Fig. 4.1b where the actual range is 1855 m and the mean is 
2050 m. 

Circles representing contours of a joint distribution of t,;t2 measurements* 
were transformed to the DR grid. The results, shown in Fig. 4.1d demonstrate the 
asymmetry of the DR measurement noise. 

The inversion bias was studied by Moose [11] for the idealized homogeneous 


case. The combined opposing effect of the inversion and the nonnegative delay are 


* Ot, = O1,;t, and tz independent. 
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Fig. 4.1. Nonlinear inversion bias. 
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examined here for the IH case. Means to compensate for these biases are discussed 


as well. 


B. THE PREFILTER EVALUATION SCHEME , 
The evaluation scheme shown schematically in Fig. 4.2 involved the following 


steps. 


evaluated case reference case 








= 


Q INVERSION DR. 






Fig. 4.2. Inversion performance evaluation. 


1. Generation of the reference and evaluated direct function tables using the 
SMART program and the instrumentation model (see Chapter Three.). Two 
separate sets of conditions are used: the “actual” (the reference case or Re- 
fease) and the “assumed” (the evaluated case or Evalcase). The conditions 
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such as SVP bottom depth grid size etc. may be different for the Refcase and 
the Evalcase. It is the effect of these differences that is being evaluated. 


2. Generation of the delay measurement noise for every point in the grid and 
adding it to the delays t,,t2, of the (reference) TDOAs. 


3. Simulation of the noisy target predicted position D,, Ry by adding depth and 
range noise to every point on the grid values Dj), #,;). This supports investi- 
gation of the search algorithm in an open loop mode. 


4. Transformation of the measured ¢,,t2 into depth and range measurement 
Dm,Rm via the inversion procedure (prefilters 5 or 6) using the evaluated 
direct function (i.e., the assumed ocean) and the noisy reference point. 


5. Comparison of the measured D,,,Rm to the reference Dj), Rij) values and 
generation of the errors 


Deer = oe oa D{t] 
iteerar = Rm 7 Rij| 


6. Repeating the above for all the points in the DR grid. Recording for each point 
the errors along with the search performance measures (number of solutions 
and search area size). 


7. Averaging the output over small local 2-D windows to provide an estimate of 
the local mean and variance and a local prefilter performance measure. 
Recall that a failed search will return the reference point as the inverted out- 

put. When noise is not added to the reference point this will produce a zero error 
which is not a precise account for the event. This will become evident from some 


of the results. 


C. EVALUATION RESULTS 
The evaluation was performed about 150 times for a variety of conditions. The 
results of eight summary cases are presented next. The cases investigate the effects 


of the following issues: 
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e Depth range grid size. 

e Linear versus bilinear interpolation. 

e TDOA bias error. 

e Acoustic medium inhomogenuity. 

e Range independent TDOA noise. 

e Range dependent TDOA noise. 

e Range dependent prediction point reference. 


e Bottom depth error. 


Each of the next eight subsections describes one evaluation. Each subsection 
describes the purpose of the simulation experiment and the expected results, the 
actual experimental results and analysis of the evaluation, and conclusions. The 


output plots include: 
e Local average of depth and range errors. 
e Local average of inversion error STD. 


e Local average of the solution count (this monitors the performance of the 


search algorithm. ) 


e Average error as function of range (range error profile). This is generated by 
taking error averages over all of the depth points for a given range. The error 
profile is normalized as marked on the plots. 


1. Effects of Grid Size 
The purpose of this Run was to examine the effect of grid size on interpo- 
lation error and optimize the grid for the intended use. The key parameters used 
are given in Table 4.1. | 


The depth and range results of the inversion with linear interpolation over 


the medium density DR Evalcase grid ($1111), are compared to the exact depth 
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TABLE 4.1 
INTERPOLATION EVALUATION GRIDS 


Grid density High Medium Low 
Casename $1155 $1111 $1144 
Initial depth [m] 160 160 160 
Depth Step size [m] 10 50 100 
Number of steps in 11 6 
Initial range [m] 500 500 1000 
Range step size [ml] 250 500 1000 
number of steps 99 20 13 


and range values of the high density DR Refcase ($1155) grid. The results are 
shown in Fig. 4.3. The maximum errors were 16 m in depth and 100 m in range. 

An interesting point was revealed by these plots. The interpolation er- 
rors depend on the degree of curvature in the direct and inverse functions. The 
curvature is very small at long ranges —s" the interpolation error inversely 
proportional to the range. This is clearly seen in Fig. 4.3. The significance of 
the interpolation errors is greater at short ranges than at long ranges since the 
other inherent MP tracking errors are reduced at short range. The grid size should 
therefore be optimized for the short range performance. 

A more esas grid (Evalcase $1144) was evaluated using the medium 
density grid as the reference (Refcase $111). The results shown in Fig. 4.4 reach a 
maximum depth error of 20 m and a maximum range error of 200 m. 

Conclusion: The grid size depends primarily on the direct function. A 


more curved and more rapidly changing function will require a finer grid to maintain 
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Fig. 4.3. A 50 x 500 m DR grid interpolation error. 
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Fig. 4.4. A 100 x 1000 m DR grid interpolation error. 
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low errors. For the nominal homogeneous case, in water depth on the order of 500 
m and ranges up to about 20 km a direct function Q grid size of 2 x 20 x 40 will 
be required. This is a very reasonable size which is easily within the capability of 
current onboard microcomputers. 

2. Bilinear Interpolation 

The purpose of this evaluation was to compare the performance of linear 
and bilinear interpolation. The key parameters used are the same as those used for 
the medium grid evaluation (See Table 4.1) with $1155 as the Refcase and $1111 
as the Evalcase. 

Error surface for the bilinear interpolation are plotted in Fig. 4.5. The 
difference between the results of the two interpolation methods is small (compare 
to Fig. 4.3). The smaller depth errors in the linear interpolation result from the 
fact that the triangular area of interpolation is smaller than that of the rectangular 
area used by the bilinear method. The bilinear interpolation provides a smoother 
interpolation but has a larger absolute error. 

Conclusion: Both linear and bilinear interpolation produce satisfactory 
results. The error of the linear method is slightly smaller. The choice of which 
method to use can be based on implementation considerations. Linear interpolation 
was used for most of this study. 

3. The Effects of Acoustic Medium Inhomogenuity 

The purpose of this evaluation run was to examine the effects of the 
medium inhomogenuity. A SVP with a positive surface gradient of steepness that 
is characteristic of many practical problems (g = 0.05 sec™*.) was used as the ref- 
erence case (A2251) and the TDOAs were converted to depth and range assuming 


a homogeneous medium (straight line propagation, Evalcase $2251). 
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Fig. 4.5. Bilinear interpolation error. 
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Depth and range inversion errors are shown in Fig. 4.6 as surfaces and as 
contour plots. Very large errors are developed especially at long ranges. Observe 
that at a range of 12 km the range error is 4 km and the depth error is 150 m. 

The trend of the error is expected and is caused primarily by the change in 
the surface bounce path. Fig. 4.7 shows a comparison of an IH and a straight line 
propagation where t, and tz for both cases are the same. The path 7; is longer 
(curved instead of straight) and includes portions with average speed of sound 
lower than in the homogeneous case. The resultant delay ¢; given by T, — Zo, is 
longer. 

When interpreted as if it were the result of straight line propagation it gives 
rise to deeper depth and shorter range estimates. 

Conclusion: The inhomogenuity of the ocean medium has a gross effect 
on MP propagation and position measurement. Compensating for this effect is 
mandatory in order to achieve meaningful depth and range measurements. 

4. Finite Delay Resolution 

This case examined the effect of finite delay resolution combined with the 
effect of the inhomogenuity. The cases used were Refcase C4261 for which the 
assumed delay resolution is 0.5 ms and Ewvalcase $4261 for which the assumed 
delay resolution is perfect and the assumed SVP is homogeneous. 

The errors are similar in nature to the errors produced by the previous 
run. The main difference is along the surface as shown on Fig. 4.8. This differer<e 
results from the fact that in the finite resolution case the surface bounce is not 
resolvable along the first depth line of the DR grid. The first depth grid point is 


located at a depth of 2m from the surface*. The delay between the direct path and 


* The observer depth is 162 m and the grid initial depth is +160 m relative to 
the observer. 
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Fig. 4.6. Effect of medium inhomogenuity (A2251). 
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Fig. 4.7. Inhomogenuity error. 
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the first surface bounce is smaller than 0.5 ms and therefore is not resolvable. As a 
result the modified TDOA vector ¢;, t2 is comprised of 7,, 73 here (instead of 7, 7). 
This corresponds to a physically different multipath configuration*. Similar effects 
were demonstrated near the ocean bottom and along the DR loci where ¢; = fg. 

Note also that the magnitude of the error over the complete DR grid is 
larger than in Fig. 4.6. This results from the steeper surface gradient of 0.105 
[sec—!] (SVP 4 in Appendix-E) used here compared to the 0.05 [sec~*] gradient 
used in Fig. 4.6. 

Conclusion: The finite delay resolution has a significant effect on the 
TDOA structure. Accounting for the specific response of the delay estimator to 
the loss of resolution is mandatory, especially for targets near the ocean boundaries. 

5. Delay Estimation Bias 

The purpose of this evaluation is to examine the effect of a constant TDOA 
offset. This was done as a scaling run and supported the independent evaluation of 
the tracker. It allowed nominal transformation of the predicted delay noise to depth 
and range measurement errors. Since this is a nominal scaling, the homogeneous 
case ($3252) is used both as the Refcase and the Evalcase. A 7 offset of 1 ms is 
used as TDOA bias. 

The resulting error is shown on Fig. 4.9 and reaches levels of up to 4 km 
in range and 80 m in depth at range of 25 km. With actual range in the vicinity 
of 12 km the range error is 1 km and the depth error is 40 m. 

Conclusion: If the tracker operates over a range of up to 12 km and 
TDOA noise of the order of 1 ms is assumed, the range error STD can be expected 
to lie between a few hundred meters and about 1 km. As was shown in Chapter 
Two, such errors can be handled by the MM target tracker. 

* See Chapter Three, Section III.B.5. 
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Fig. 4.9. TDOA constant bias error (7; + lms). 
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6. Effects of Range Independent Noise 

The effect of a stationary (range independent) random noise on both the 
variance and the mean of the DR measurement was examined. The Refcase and 
Evalcase were both taken as the homogeneous 54261. A stationary, range indepen- 
dent ((p = 0), see Eq. (3.3)) TDOA noise with STD of 1 ms for both ¢; and t2 was 
used. The local averaging window for the plots was 3 x 7 in grid step units. 

Surface plots of locally averaged STD of the depth and range errors are 
shown on Fig. 4.10a. The depth and range errors were averaged over depth to 
provide the range error profile which is shown on Fig. 4.10. The profiles are 
normalized to 100 m in depth and 1 km in range. 

The maximum error is consistent with the earlier scaling results performed 
using the constant TDOA offset (Fig. 4.9). Contours of locally averaged number 
of solutions (Fig. 4.10d) indicate frequent failure of the inversion near the deep 
end of the grid*. This is caused by the TDOA values produced from the noisy 
measurements. The noise drives the TDOA measurements out of the interval of 
TDOA values covered by the precomputed direct function. Increasing the area 
covered by the DR grid will solve part of the problem and reduce the frequency of 
failures. 

The results of a similar case with TDOA error STD of 5 ms are shown in 
Fig. 4.11. This figure shows larger errors and failure of the algorithm at a shorter 
range (7 km in Fig. 4.11d). Note that the term failure indicates only that the 
effective sampling rate is reduced by half at ranges over 7 km (recall the comment 


at the end of Section B). More frequent failure is the reason for the superficially 


* The frequency of failure is 50 % along the 0 solution contour line since the 
line separates the regions between 0 and 1 average solutions. The contour line does 
not imply that there are no solutions at ranges longer than those indicated by the 
contour. 
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Fig. 4.10. Range independent TDOA and DR noise (1 ms). 
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lower error values at the long range. The effect of the reduced sampling rate on 
the overall system is investigated in Chapter Five. 

Conclusions: Errors with STD on the order of 30 m in depth and 0.5 km 
in range are expected at ranges around 10 km as a result of a range independent 
TDOA noise with STD of 1 ms. For TDOA noise with STD of 5 ms similar depth 
and range errors are expected at a range of 5 km. As shown in Chapter Two these 
errors can be handled by the target tracker. The inversion algorithm does fail at 
long ranges if the noisy delay measurements fall outside of the interval covered by 
the precomputed direct function. 

7. Effects of Range Dependent Noise 

The response of the inversion of TDOA which is contaminated by range 
dependent noise was examined here. Range dependent noise was added to the 
target predicted point as well. Case $4261 was used both as the Refcase and the 
Evalcase. Range dependent noise (Eq. (2.96) NT =1, P = 1), with t),t2 STD of 
0.2 ms at initial range of 500m, was used. The TDOA noise STD varied linearly 
with ranges between 0.2 ms to 4.8 ms over the range interval of 0.5 to 12.5 km. 
This noise level resembles the one produced by the more elaborate noise model 
(NT=2) for a SNR of 60dB and p = 2 which is a fairly typical case. Predicted 
position DR noise was varied from 0.4 to 250 m in depth and from 10 to 6.25 m in 
range proportional to the second power of range (ps, = 2). 

The resulting error profile is shown on Fig. 4.12 along with the contour 
line for the number of solutions. Note the sharp range dependence of the inversion 
error STD on range. Also note that for the practical case used here a sharp “knee” 
develops around 12 km (Fig. 4.12c). Beyond this range the error grows very 
rapidly. This close to second order dependence of the DR measurement noise on 


range 1s what motivated setting the parameter p,; to 2. The maximum area size 
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Fig. 4.11. Range independent TDOA noise (5 ms). 
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was set to 20 x 25 grid points. This explains why the addition of the DR noise did 
not degrade the performance significantly. Performance degradation as a result of 
smaller maximum area size is demonstrated in Chapter Five. 

Note that the normalized range STD in Fig. 4.12b at a range of 9 km 1s 
approximately 0.3. This implies a range error STD of 3.6 km (considering the 50% 
failure rate and the 6 km normalizing factor). This value will be useful in Chapter 
Five in setting the maximum area size. 

A second interesting observation is the failure along the first range grid 
line. This is caused by the noise which drives the “measured” TDOA outside of 
the interval covered by the precomputed direct function. Recall however that the 0 
contour line in the number of solutions plot indicates only that the average failure 
rate is larger than 50%. Similar phenomena were noticed in most simulations along 
the edges of the DR grid. 

Conclusions: The inversion varies with range even for a range indepen- 
dent TDOA noise. The error becomes even more range dependent when the TDOA 
noise itself is range dependent. Beyond 15 km very large DR measurements develop 
even for moderate initial range delay estimation noises. 

8. Effect of Bottom Depth Error 

The purpose of this last inversion evaluation was to quantitatively examine 
the effects of an erroneous bottom depth on the target depth and range errors. 
Refcase $1256 for which the bottom depth is 513 m was compared to Evalcase 
91556 for which the bottom depth is 550 m. A homogeneous ocean, perfect delay 
resolution and no noise were used here in order to isolate the effects of the erroneous 


bottom depth. The inversion was done, however, using the numerical algorithm 


(prefilter-5). 
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Fig. 4.12. Range dependent TDOA noise. 
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The effect of this relatively large bottom depth error (37 m which is 7.2% 
of the total ocean depth) is shown on Fig. 4.13 for short and long ranges. While 
the resulting range error reached a level of 3 km at 18km, it was limited to 1.5 km 
for ranges smaller then 10 km. 

The trend of the error surfaces is explained as follows. The actual (Refcase) 
ocean produces a bottom bounce TDOA which is shorter then the one produced for 
the same DR point by the assumed (Evalcase) ocean. This smaller t2 is interpreted 
by the inversion ( which is done using the Evalcase) as the TDOA of a deeper target 
at a longer range. The region of no inversion solution at ranges in excess of 18 km 
resulted from the Refcase tg TDOA being smaller then the minimum f?2 in the 
precomputed Evalcase Table. 

Conclusions: Incorrect bottom depth assumptions can give rise to large 
errors in MP position measurements, primarily at long range. The trend of the 
error is consistent with other MP tracking error sources in that it increases with 
increasing range. This error tends to limit practical tracking in shallow water to 


ranges below 20km. 
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SUMMARY OF PREFILTER EVALUATION 


The results of these experiments can be summarized as follows. 


. The effect of the inhomogenuity on the MP position measurement is severe 
and was clearly demonstrated. 


. The depth and range measurement errors are strongly dependent on range 
(proportional to 2"¢ to 4** power of range depending on SNR and propagation 
loss). 


. Errors with STD of up to 40 m in depth and 1 km in range can be expected 
for ranges of about 7 km and typical targets with SNRo = 60 dB if spherical 
spreading is assumed (p = 2). Above 12 km the errors increase sharply with 
range for the typical shallow water case studied here. 


. The inversion algorithm fails under very noisy conditions and around the edges 
of the DR grid. The failures in these cases are legitimate since the TDOA input 
values are not representative of points in the DR grid. 


. Errors will result 1f account is not taken for the loss of delay resolution. In 
particular one needs to model the response that is expected from the specific 
type of delay estimator. 


. The algorithm is sensitive to DR grid step size only at short ranges. In general 
the depth and range step sizes should be tailored to the fluctuation of the 
specific direct function used and the accuracy desired at short range. Grids of 
10 to 50 m in depth and 250 to 500 m in range were used in most of our work. 


. The MP tracking is sensitive to bottom depth errors primarily at long ranges. 
At ranges below 10 km a low percentage error in bottom depth would produce 
about twice that percentage error in range. 
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V. OVERALL PERFORMANCE EVALUATION 


A. INTRODUCTION 

In this chapter the overall performance of the MP tracking system is evalu- 
ated. A total of seven simulation runs which are divided into into three groups 
are presented. The runs are numbered 6 through 12 following the sequence of the 
runs begun in Chapter Two. Runs 6 and 7, represented in section B, demonstrate 
the functional performance and limitations of the integrated prefilter and tracker. 
Runs 8 through 10, described in Section C, were designed to investigate the effects 
of medium inhomogenuity and the ability of the new prefilter to compensate for 
it. Runs 11 and 12, covered in Section D, examine the sensitivity of the IH com- 
pensation to the accuracy of SVP measurement. Conclusions are summarized in 
Section E. 

The results are presented in a format similar to the one used for the previous 
runs and case evaluations in Chapters Two and Four. In each section the goal 
of the runs and the key parameters used are described first. Then the analysis 
results are presented. The detailed parameter values for all of the runs is given in 


Appendix E. 


B. INTEGRATED SYSTEM FUNCTIONAL PERFORMANCE 

The overall performance of the prefilter integrated in a closed loop with the 
target tracker was demonstrated using a maneuvering target. The target initially 
has an outgoing speed (range rate) of 6 m/sec and it makes a sharp maneuver at a 
range of 9 km turning back towards the receiver at the same speed. Realistic noise 


(noise type 2, SNR = 50 dB and p = 2) is used to investigate the overall 
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capability of the system through the maneuver. A constant speed of sound is used 
both as a reference and evaluation medium (case $1111). This is done in order 
to concentrate on the functional performance of the integrated system when the 
measurements are noisy. 

The key issues investigated here relate to the inversion method. These issues 
are the ability of the local search algorithm (LSA) to find a valid surrounding 
triangle and the effects that a reduced effective sampling rate has on the target 
tracker (recall the discussion of inversion failure in Chapter Three). 

Results of Run 6 are shown in Fig. 5.1. The maximum search area size was 
+ 200 m in depth and +3.75 km in range around the predicted target position 
(a rectangle of 8 x 15 DR grid units for case $1111). Fig. 5.1 shows that the 
overall tracking performance is good. The range tracking error reduces to around 
300 m at the 14** minute after initialization is completed. The realistic TDOA 
noise simulated here increases rapidly with range. This is clearly evident in Fig. 
5.1 (note the two resulting range measurements errors around the 13** minute). 

The system tracks well through the maneuver with an effective sample rate as 
low as 40%. This is indicated by the average number of solutions which is close 
to 0.4 in Fig. 5.1d. Note that the prefilter uses the maximum search area size 
during the 17" minute as indicated by the normalized size in Fig. 5.1d (actual 
searched area / maximum area size). Some ambiguity is indicated at shorter ranges 
where the solution count is larger than 1. This ambiguity affects mostly the depth 
tracking as can be seen in Fig. 5.1b. This effect is expected when one reviews 
the TDOA contours of case $1111 which is used here (see Fig. 3.6). The bearing 
channel tracking and the resultant XY plot of Run 6 are shown for completeness 


in Fig. 5.2. 
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Fig. 5.1. Run-6: IH tracking of a maneuvering target. 
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Fig. 5.2. Run-6: Bearing and XY plot. 


The same scenario was used in Run 7 with a smaller maximum search area size 
of +125 m in depth and £2.25 km in range (a rectangle of 5X9 in the $1111 grid). 
The results, shown in Fig. 5.3, clearly indicate the failure of the system to track the 
target through the maneuver. The system loses track when the range tracking error 
reaches 2 km (in the 20°" minute). At this time and range most of the measurements 
are outside of the maximum search area around the predicted position. This leads 
to a constant use of the maximum search area size and eventually to a complete 
failure (over 8070 of the iterations). 

The range measurement error STD for Runs 6 and 7, was estimated using the 
inversion performance evaluation scheme described in Chapter Three. The typical 


STD for similar TDOA noise at around 9 km was 3.6 km*. The conclusion is 


* Recall the discussion in Chapter Four, Section C.7 
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Fig. 5.3. Run-7: Failure due to reduced search area size. . 
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therefore that the search area size should be set to twice the standard deviation of 


the measurement noise, ((2¢4) x (2¢,)). 


C. TRACKING IN AN INHOMOGENEOUS MEDIUM 

The next three Runs demonstrate the ability of the new MP tracking scheme 
to compensate for the effects of ocean inhomogenuity. A reference case (Refcase) 
was used to generate the TDOAs (the direct function table) used for simulation of 
the “actual” medium. The precomputed direct function table generated using the 
evaluated case (Evalcase) was used in the prefilter. The Evalcase represented the 
assumed ocean conditions. 

Run 8 demonstrates the effects of an IH medium (Refcase C2251) when straight 
line propagation is assumed and used in the prefilter (Evalcase $1111) . The results 
are shown in Fig. 5.4. The large range tracking error of 3.3 km at a range of 10 km 
clearly indicates the need for compensation for the medium inhomogenuity. A 
relatively strong target signal (SNR» 70 dB) was used in order to isolate the IH 
medium effects from the effects of the noise. 

Run 9 examines the same scenario with a weaker target signal (SNRo = 50 
dB) and its results are shown in Fig. 5.5. The TDOA measurements are noisy now 
and give rise to a negative range bias which adds to the effect of the IH ocean. The 
total resulting estimation error for this case is 4 km at the 10 km range. 

The use of the new prefilter with the correctly assumed IH medium (Evalcase 
C2251) is demonstrated in Run 10 and the results are shown in Fig 5.6. The 
removal of the error associated with the IH medium is clearly evident. The error 
in the depth channel after the depth maneuver at the 8'* minute is caused by 
the fact that the C2251 grid was limited to a depth of +60m (due to the current 


limitation of the SMART eigenray program discussed in Appendix F). This error 
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is not present at shallower depths away from the grid edge or if the noise is smaller 
as demonstrated by the near perfect tracking in Fig. 5.6 where SNRog of 70 dB 
was used. This error can be removed completely if a full grid 1s employed. The 
average number of solutions reached 1.8 in some cases since the corresponding 
TDOA contour lines intersect in four closely located DR positions (see Fig. 3.17 
which is taken from typical contours of the same case C2251). The effects of the 
increased TDOA measurement noise are evident at the long range toward the end 
of the run. 

The overall tracking capability in an IH medium demonstrated by this run is 
also shown in Fig. 5.6 as an XY plot. The estimated track (heavy line) approaches 
the actual track (marked with dots) whereas the measurements are very noisy (the 


fluctuating line). 


183 


RANGE [KN} 


RANGE [M] 
—3000 


10 


—1000 


—2000 


—4000 


-5000 


RANGE:ACTUAL *, MEASURED + , ESTIMATED V 


12 
TIME [MIN] 


RANGE ERROR:MEASURED + ESTIMATED V 


SSA 
\ \ 


nee 


4 8 
TIME [MIN] 


16 


DEPTH[ 10M] 


SIZE+MAXSIZE +, NUM OF SOLUTION ¥ 


DEPTH:ACTUAL *, MEASURED + , ESTIMATED ¥ 


15 


10 





= 


=o 


12 


0 ;% 8 
TIME [MIN] 


16 


PREFILTER PERFORMANCE 


1.0 


0.8 


0.6 


0 4 4 
TIME [NJ 


12 


Fig. 5.5. Run-9: Effect of IH and TDOA noise. 
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Fig. 5.6. Run-10: Correction for medium inhomogenuity. 
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D. EFFECTS OF ERRONEOUS SVP 

The sensitivity of the IH MP inversion to errors resulting from an assumed 
SVP which is different from the true SVP is examined in this section. A SVP with 
positive surface gradient of 0.105 sec~' was used as the actual medium in Run 
11 (Refcase A4261). A similar but slightly different SVP with a surface gradient 
of 0.100 sec~! was used in the prefilter (Evalcase A3253). The results are shown 
in Fig. 5.8. The range error is quite small for ranges less than 7 km. At ranges 
greater than 8 km the TDOA resulting Seni the true SVP using the Refcase are 
not covered by the direct function used by the prefilter. This again is not a limit 
of the basic method but rather a result of the limitation of the current SMART 
model which forced the grid selection (see Appendix F). For the region evaluated 
here one can conclude that the method is insensitive to small errors in SVP. This 
trend is expected to extend to longer ranges. 

When large SVP errors are present, the resultant tracking error is also large. 
This is demonstrated by Run 12 where a gradient of 0.105 sec~! was used in the 
Refcase and a gradient of 0.05 sec™’ was used as an Evalcase (C4261 and C2251 
respectively). The results, shown on Fig. 5.9, indicate large tracking errors that 
resemble the effect of the uncompensated IH medium demonstrated in Run 8. This 
is expected since the SVP error in the assumed SVP in Run 12 is of the same order 


of magnitude as the total uncompensated inhomogenuity in Run 8. 
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Fig. 5.7. Run 10a: Correction for medium inhomogenuity, low noise. 
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SUMMARY OF SIMULATION RESULTS 


The results of the complete system evaluation can be summarized as follows: 


. IH propagation introduces significant errors in MP tracking. The effect must 
be compensated for in order to achieve accurate tracking. 


. The IH prefilter proposed here is functionally sound and can provide compen- 
sation for the effects of the IH medium. 


. The DR grid used to compute the direct function should be chosen with care. 
It should cover all possible target positions, with an additional margin to 
enable handling of measurement noise. Too wide a grid, primarily along the 
depth axis, adds ambiguous solutions and should be avoided. 


. As a general rule, the maximum search area size should be made larger than 
twice the standard deviation of the depth and range measurement noise, 1.e., 
area size > (204) x (20,). 


. The maneuvering target tracker can handle reduced effective sampling rates 
resulting from inversion failures. The original sampling period was selected 
here as 5 sec to match target dynamics and to provide a reasonable TDOA 
observation time. Effective sampling rates which at times are reduced to as 
low as 40% of the original rate can still support reasonable tracking. 


. The inversion is insensitive to small errors between the assumed SVP and the 
actual one. (5-10% in the sound speed gradient) 


. The overall performance of the MP tracking degrades rapidly with increasing 
range. This results from the vanishing TDOA at longer range and the rapid 
increase in TDOA noise due to the attenuated target signal. 


. As aresult of the growing noise at longer ranges, measurement biases develop 
which tend to underestimate the range. Eventually the MP tracking filter fails 
completely and losses track of the target. Reliable tracking can be expected 
up to ranges of about 15 km. 
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VI. CONCLUSION 


This chapter concludes the thesis. It includes a brief review of the problem; 
a short description of the solutions and a summary of the results. It ends with a 
suggestion for further extension of the work. 

Current acoustic broadband MP tracking techniques, which have evolved over 
the last ten years, are based on straight line propagation. In reality the ocean 
medium is inhomogeneous and the sound propagates along curved ray paths. MP 
tracking therefore encounters large errors if the IH propagation is not accounted 
for. Since there is no closed form relation between the TDOA and source depth and 
range (DR) for the inhomogeneous case, the compensation for the inhomogeneity 
is computationaly complex and until now had not been done. 

Previous MP tracking algorithms have also assumed that the TDOA are al- 
ways resolvable and associable with the corresponding path. In practical tracking 
situations the TDOA are not always resolvable due to finite bandwidth of the re- 
ceiver and association of TDOA with the corresponding path is very difficult. The 
impact of these assumptions on MP tracking is significant. 

In addition TDOA measurements are noisy. The SNR depends on the strength 
of the received signal which in turn depends on range due to the propagation loss. 
A realistic overall evaluation of a MP tracking system requires a more accurate . 
model for the noise than those previously used, as well as a simulation of the IH 
medium effects. 

Current MP tracking algorithms employ the multiple model (MM) estimation 
technique. This technique produces tracking biases since the multiple models are 


in most cases not centered around a true model of the target. 


Toa 


This work focused on developing a model that accounts for the IH propagation 
and the effects of realistic TDOA processing. A key contribution of this work is the 
development of an algorithm to compute depth and range from TDOAs produced 
by MP in an IH medium. In our approach the TDOAs are assumed to be estimated 
by realistic instrumentation and thus have limited resolution and are not assumed 
to be associable with their corresponding path. 

The performance of the model is evaluated using a accurate range dependent 
model for the noise. The depth range calenlatilen is performed with the help of an 
elgenray program exercised over a DR grid. The realistic TDOAs are numerically 
computed offine and stored in a direct function table. The measured TDOAs 
are then used online to produce the corresponding DR using linear interpolation. 
The method accounts for any multipath structure and is not dependent on the 
association of the TDOA with the multipath. 

The interpolation is performed as follows. The DR and TDOA values of three 
points in the direct function table that have TDOA values close to the measured 
set, are used in the interpolation. The triangle of the three closest points in the 
direct function table is located by means of a unique search algorithm called the 
local search. The search is limited to a small area centered around the predicted 
target position in order to improve the efficiency and reduce the ambiguity. 

A technique for centering the MM hypotheses bank around the estimated 
target maneuver command to reduce the tracking bias was also developed. The 
entire model is shifted in a manner that requires no computation. It 1s based on an 
analytically derived linear relation between the command hypothesis bank and the 
states of the corresponding models. A second order smoother instead of the first 


order one previously used eliminates the lag error encountered in earlier tracking 
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designs. An evaluation of the performance of the new MP tracking scheme in 


shallow waters using a realistic model for the noise shows very good results. 


A comprehensive set of test cases were run to evaluate the components of 


the MP tracker both individually and as a total system. The main results of the 


simulation show that: 


i. 


Significant errors are introduced in MP tracking by the ocean IH acoustic 
medium. 


j 


. The new MP inversion method and the tracking technique can compensate for 


the IH medium effects. 


. The effects of limited delay resolution and lack of TDOA and path associability 


are also accounted for by the new method. 


. The depth measurement is more affected by the ambiguity resulting from the 


lack of TDOA path associability than the range measurement. 


. Realistic noise limits the performance of shallow water MP tracking to short 


ranges (say less than 15 km). The noise initially introduces estimation biases 
and eventually leads to a complete loss of track. 


. At short ranges the method is insensitive to small errors in SVP. The method 


is more sensitive to errors in ocean bottom depth. Large errors in SVP or 
bottom depth introduce large MP measurement errors. 


. Recentering the MM around the estimated maneuver command and the use 


of the second order smoother significantly reduces tracker estimation bias. 


. The coordinate decoupling technique currently used in MP tracking introduces 


tracking errors for nonmaneuvering targets around the CPA at short ranges. 


. A second order target model is sufficient for maneuvering target tracking if 


the Kalman gain is properly trimmed. 


There are several directions in which this work can be extended. However the 


following specific direction is a natural extension of the ideas and offers potential 


improvement of the results derived here. 
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Use of more than the first two TDOAs has the potential of reducing the depth 
and range estimation noise. More importantly use of additional TDOAs will reduce 
the ambiguity in target position. If, in addition, information on the reflections is 
available* and the signal spectrum is known then the complete ACF of the received 
signal can be constructed for every point in the grid. A classification method could 
then be set up to measure the closeness of the measured ACF to those constructed 
for the points in the grid. The conditional probabilities of the target position for 
any of the DR grid points can be generated by the MM tracker and incorporated 
into the classifier. Interpolation between the closer DR points can then be used to 
compute the target position. 

In summary, this research investigated the influence of realistic conditions on 
the multipath tracking of maneuvering targets. The inhomogeneity of the acoustic 
ocean medium; the limitations of TDOA estimation; and the difficulty in associa- 
tion of TDOA with multipath all have a significant effect on the tracking. A new 
model was devised to counter these effects. A very thorough evaluation of the model 
was conducted with a special emphasis on the reality of the simulated condition. 
The results indicate that in general the MP tracking is limited to short ranges. 
Accurate tracking can only be achieved if the distorting effects of the medium and 


the delay estimation are compensated for as was successfully done here. 


* Such as when operating in a well known and previously measured area 
of the ocean 
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APPENDIX A 
TRANSITION MATRICES 


A separate model of the form 
Xn = PXn-1 +TUn_-1 + YWa-1 (A.1) 


is used to describe target motion along the X,Y or range R axis. The components 


of the matrices appearing in this equation are given below. 
$i; =1 
$21 = $31 = $32 = 0 
$12 = (1— 7°") Ja 
$22 = ep= 


oi. = [1 + (aye % — Ger) (a - Aw )| — 








1 

- Wl ae \ 
$23 = (e e€ ) i _ 
$33 = owt 
P, = (aF -1+e7°7) = 
te —l—en 
rs =0 

= Ow _ ,a~al eo _— way I 1 
v= {T+|% (1 e ) “a (1 e =} 

1 

—_ 2. pw Owl \o— _ pwal 

VY, = a, (1 e ) ol e ——- 


| 
YY, = (1 -- ence.) ray 
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These parameters are easily derived from Eq. 2.3.3 in Ref. 6 by the change of 


variables 
1 
Wane —_ (A.3) 
AQw 
and 
Gai G. (A.4) 


The advantage of the scaling used here is that both U and W are now in [m/sec], for 
example a constant process noise W or command U of 1 [m/sec] in the X direction 
would each yield z = 1 [m/sec] at the steady state. 


For the cross range channel the equation is 


B= $5: Ba-1 +T5Us + Wown-1 (A.5) 
where the components of the matrices are expressed 


in terms of the previous models components as 


$11 = $11 

$21 = $531 = $o32 = 0 
Pri2 = Piz 

P22 = rr 

$o31 = $31/R 

$532 = $32/R 

$433 = $33 

Tee Di/ 2 

Ty. =T2/R 
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I',3 =0 


V,,=V,/R 
V,. = W2/R 
v,; = V./R 


These results derive from Eq. 2.3.28 in Ref. 6 with the change of variables Eq. 
(A.2) - (A.4). — 
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APPENDIX B 
CONDITIONAL MEAN ESTIMATE 
A constant U is to be estimated from the observation z = U +n where 
n ~ N(0, 07). The optimal mean square estimate is the conditional mean U = E{U|z}. 
Two equiprobable hypotheses U,,U2 are assumed, which span the range of the ex- 


pected U. The estimate is thus given by 


U = E{U|z} = su, - P(U;|z) = yu, aa (B.1) 


where f(z) is the density function. 
With P(U,) = P(U2) = 1/2 and the given normal distribution of the noise 
the estimate U is given by: 


Pd race ele U,)?/20? cee ve e—(2—-U2)* /20° 








B.2 
Te ee = Cariece ale Faze an 2—U2)-/ 20" ( ) 
which after some algebra is: 
oe xi Xx 
_ 1 + e—[—-22(U2 —U,)—(U2-U;)(U2-U )] /2e0? a 1 + e—e{2z(U2 —U,)—(U2 —U,)(U2+U; )] : 
(B.3) 
If we let 
U, =U, —-U,;U, = (U, + U,)/2; (B.4) 
and 
Us 
a = a2 (2 = oe) (B.5) 
then we obtain 
F Fe Ws OP . U,e7%/? (Eon fe eT ue | ite Uz e%!? (enor 4 e— 7 
l+e*” t1+e-¢ (eo/2 4. e-a/2)" 
(B.6) 
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‘ which after substitution of Eq. (B.4) becomes 


(U. —U,/2)e—*/? + (Ue + Us /2)e%/? _y AD B23 gle _y Os mee 
(e%/2 + e—@/2) _  ) 64a/2 4+ e—a/2 B bs) 9° 


(B.7) 


U = 


After substitution of Eq. (B.5) the result is 


. U, Us 
U =U. + = tanh oe: (z 
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APPENDIX C 
COMMAND UPDATE TRANSIENT 
In this appendix the transient in command estimate as a result from a com- 


mand bank update is computed. The transient is defined as 


U pee — U.+ = Vine (Oa) 
and is given by 
Urania = UI,- ° W,- <a UI,.4+ ; W y+. (CZ 


The transient is the result of the command bank update given by 
UI,+ = UI,- +N- AU. (C.3) 
Examining Fig.2.8 one realizes that for: =2 ton 
(Ui )e- « (Wide- = (Ui-1)at + (Wi-1)a+- (C.4) 


Applying the above to Eq. (C.2) cancelling terms and recalling that (W; ),.- = (Ww )xz+ 


results in 
Utransient = (Un) a+ «(Wn )et — (Ui )a- -(Wi)e-- (C.5) 
But since 
(Un e+ =(U1),- +N -AU (C.6) 
then 
Vicenwent= N i - GPip sa. (C.7) 


For an update in the opposite direction one gets with the same arguments 


— =—-N-AU.- (Wn )x- (C.8) 
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APPENDIX D 
DELAY ESTIMATION ERRORS 
The transformation of the acoustic noise to a biased range dependent delay 
noise by the time delay estimator is considered here. This noise modeling was 
done in order to provide a more accurate simulation than previous MP tracking 


simulations. 


A. VARIANCE DEPENDENCE ON SNR 

The minimum error variance in estimation of a variable from a measurement 
with zero mean random noise is given by the Cramer Rao lower bound (CRLB). 
The bound has been applied to TDOA estimation. Ianniello [Ref 14] gives a cur- 
rent and relevant review of time delay estimation error. The performance of an 
autocorrelator time delay estimator is evaluated there both analytically and exper- 
meettaty. 

Iannoello shows that if the z(t) signal and the noise n(t) both have uniform 
spectral density in the band 0 - B Hz and are zero outside of that range, then the 


variance of the autocorrelation estimate of r from the single multipath signal 
y(t) = z(t) + az(t — 7) + n(t) (coe) 


where a is the relative attenuation (gain) of the path is given by 


: (1 +a?) + a | S2+2(1+a?)Sg-Ng + N2 
———————— Se ee 


aan (D.2) 





(D.3) 
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where S and N are the spectral densities of the signal and the noise in the 
baseband 0 — B and T is the observation time. For the special case of a = 1 the 
variance reduces to 


3 —l —2 


which has three dominant regions depending on the signal to noise ratio 
SNR = S/N. (D.5) 


For SNR >> 1 the variance approaches a constant indicating the performance limit 


of this type of estimator. 
| ae (D.6) 
For SNR << 1 the error variance is inversely proportional to the second squared 
SNR. In between these two regions there is a transition from constant to dependence 
on SNR7?. 
In his paper Ianniello indicates close agreement of Eq. (D.2) both with sim- 
ulation results and with the CRLB for the typical case of SNR << 1 encountered 


In practice. 


B. VARIANCE DEPENDENCE ON RANGE 
For a stationary acoustic noise with bandlimited spectral density NB and signal 
intensity given by Eq. (3.3) the SNR at the input of the delay estimator becomes 


range dependent and is given by 
oN (y= SN, -7F (D.7) 
where SNR, is the SNR at range R = 1 [ml]. 
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Two models were used here to represent the range dependence of the noise. 
The first* was a simple dependence on some power p of range within a constant of 


proportionality a9 selected to match realistic noise at ranges 5 < R < 15 km. 
C= d0° RP (D 8) 


This model was used to develop the general trend of the inversion and the tracker 
when the noise is range dependent. The second noise model* was based on the 
exact Eq. (D.2) with SNR being the function of range given in Eq. (D.6). Some 
examples of the STD of TDOA as function of range for typical values of SNR arc 
shown in Fig. D.1 and D.2. 


C. DELAY ESTIMATION BIAS 

Two of the sources contributing to MP bias errors are associated with the 
delay estimation. These are the nonnegative character of the delay estimate and 
the effect of close ACF peaks. These are discussed next. 


The symmetry of the autocorrelation function, namely 
ACF(r) = ACF(-r) (D.9) 


prevents the discrimination between positive and a negative time delay (a lagging 
versus a leading replica). Since in the multipath situation, there really is a first 
arrival followed by lagging replicas, the use of the positive part of the time lag axis is 
conceptually justified. This however implies that the error in time delay estimation 
is not a simple gaussian, but rather reflects around the r = 0. A distribution which 


is equal to 0 for negative time delays, is used here to describe the estimate 7. An 


* Noise type (NT) one. 
* NT two 
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Fig. D.1. Range dependent TDOA noise type 2 (low). 
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Fig. D.2. Range dependent TDOA noise type 2 (high). 
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estimate r. An estimate with this distribution 1s generated by taking the absolute 
value of the sum of the time delay and a normal zero mean noise. 

When the actual time delay 1s large relative to the estimation error variance, 
the effect of reflection around the 0 lag is marginal. When the actual delay is small 
compared to the error variance the resulting error bias has a stronger effect. 


Figure D.3 shows three distributions of the process 
t= |t+n| (D.10) 


for t taking the value of 1.80, 0.850,0.220 where n ~ N(0, 07). It is evident that 
when ¢ is smaller than 0.50 it will be significantly overestimated which will cause 


an underestimated range. 
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Fig. D.3. Nonnegative bias. 
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D. CLOSE ACF PEAKS 

When the actual delays are short compared to the signal correlation time (1/B) 
the presence of an adjacent peak will shift the location of the estimated peak. This 
shift can be either positive or negative depending on the specific set of MP and 
resulting ACF. Since it is analytically predictable it can conceptually be corrected 
for as is done in the maximum likelihood estimator [14]. This error was therefore 


not simulated here. 


207 


APPENDIX E 


DETAILED DATA DEFINITION AND SELECTION 


This appendix describes in detail the parameters used in the various simula- 


tions and evaluations. Run data parameters are described in Section A. The values 


selected for the 10 simulation Runs are listed in Section B. Section C describes the 


case code name. 


A. SIMULATICIN RUN DATA 


Target and Scenario Data 


Ad - 


AT - 


Pln - 


Medium 


Case - 


Mdn - 


Pin - 


Pp = 


[sec] The target’s depth time constant (1/aq in Eq. (2.5). 
[sec] The target’s horizontal time constant (1/a@ in Eq. (2.2). 
The pilot number 1, 2 or 3 per Chapter Two, Section E. 
and Prefilter 


A selection of parameters defining the ocean medium 
and its acoustic modeling. 


Medium number 1! through 4. 
Prefilter number 1 through 6. 


[number] Power of the range dependence of the 

delay estimation variance (See Chapter Two, 

Section D.4, Chapter Three, Section A, and Appendix D.) 
[number] Power of range dependence of bearing 
estimation noise, (and of DR noise when it is 

used with medium and prefilter number 1, Eq. (2.101)). 
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Two cases were used for each Run. 


Refcase - The reference is a set of parameters considered as 
the “actual” conditions used by the medium 
program to convert the actual DR position to TDOA. 


Evalcase - The evaluated ocean conditions and used by 


the prefilter in the inversion of TDOA to DR. 


NT - heise type 1 = the simplified model Eq. (2.96). 
2 = the improved model, Appendix D. 


Or - [ms] Standard deviation of TDOA 
t; (or 7,) estimate at range of 500 m. This parameter 
is was used only by noise model 
type 1, but is computed for noise model 
type 2 for reference as well. 


O42 - Same as above but for t2(72). 


Od - |{m] Depth error standard deviation at range of 
500 m, used with medium 1 and prefilter 1. 


oF - |[m] Range error standard deviation at range of 
500 m, used with medium 1 and prefilter 1. 


Oo - [deg] Bearing standard deviation at range of 500 m. 
Bw - [Hz] Receiver bandwidth. 

Atn - [dB] A vector of surface and bottom bounce losses. 
SNRo - [dB] Signal to noise ratio for target at range of 1 m. 


measured at the input to the delay estimator. 


MM Parameters 


a - [sec] System sampling time interval. 
Ord - MM target model order, 2 or 3. 
Aq - [sec] Command noise coloring time constant (1/a, Eq. (2.7)). 
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Genre - [m/sec] Lowest speed command in the hypothesis bank. 
ea - [m/sec] Highest speed command in the hypothesis bank. 
N - Number of filters in the MM. 
er - Command transition probability 0 to 1. 
Anl - W averaging constant, (iterations to reach 637%) Eq. (2.23b). 
An3 - Wop, Uop averaging constant, (iterations to reach 63%) 
Eq. (2.26) and (2.27) 
An4 - Xop averaging constant, (iterations to reach 63%) 
Eq. (2.28) 
An - U, averaging constant, (iterations to reach 63%) 
Eq. (2.88) 
B. Runs 


The specific values of the parameters for each run are given below: 


Run 


- Jl 


Target and scenario data 


Ad 
Ar 
Pln 


Mdn 
Pfn 


Pb 
Od 
Or 
Ob 


MM 


7 


ord 


- 35 sec. 
- 40 sec. 


parameters 
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Aw 
oat 
Citi s 
N 


Pr 

Anl 
An3 
An4 
Ang 


Run 


Pln 


Mdn 
Pin 
Refcase 
Evalcase 
NT 


Pb 
Ort1 


Ot2 
Ob 
Bw 
Atn 
SNR 


Run 


Ot1 
O12 


SNR 


Run 


Ot] 


25 [sec] 
—30 [m/sec] 
+30 [m/sec] 
11 

0.9 

10 

4 

+) 

100 


2 (All parameters are the same as for Run 1 
except the following:) 


3 


3 

3 

sy LIL 
Saale 

y 

1 

1 

0.005 [m] 
0.004 [m] 
0.08 [deg] 
2000 [Hz] 
5.5 [dB] 
50 [dB] 


3 (all the parameters are the same as for Run 2 
except the following:) 


2 

0.05 [ms] 
0.05 [ms] 
70 [dB] 


4 (All parameters are the same as Run 2 except 
for the following: ) 


z 
0.06 [ms| 


yAal 


O12 


SNR 


Run 


Ot} 
O12 


Atn 


Run 


Pln 


Mdn 
Pin 


Pp 

Pb 

Ctl 

Ot2 

max area size 
Atn 

SNR 


Run 


max area size 


Run 


Refcase 
Evalcase 
NT 

P 

Pb 

Ot) 

O12 


0.08 [ms] 
60 [dB] 


5 (All parameters are the same as for Run 2 
except the following:) 


2 

0.018 [ms] 
0.03 [ms] 
5.10 [dB] 


6 (All parameters are the same as Run 2 
except for the following:) 


1 


to OT 


1 

0.018 [ms] 
0.018 {ms] 
8x15 

5.5 [dB] 
50 [dB] 


7 (All parameters are the same as Run 6 
except for the following:) 


ox9 


8 (All parameters are the same as Run 6 
except for the following:) 


C2251 
51111 

Z 

Z 

2 

0.005 [ms] 
0.005 [ms] 
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Atn 
SNR 


Run 
Atn 
SNR 
Run 
Refcase 
Evalcase 
Run 
Refcase 
Evalcase 
Run 
Refcase 
Evalcase 


O11 
O12 


Atn 
SNR 


55 [dB] 
70 {dB} 


9 (All parameters are the same as in Run 8 
except for the following:) 


5.10 [dB] 
50 [dB] 


10 (All parameters are the same as in Run 9 
except for the following:) 


C2251 
C2251 


11 (All parameters are the same as in Run 10 
except for the following: ) 


A426 
A3253 


12 (All parameters are the same as in Run 11 
except for the following:) 


C4261 
C2251 
0.005 [ms] 
0.005 [ms] 
5.5 [dB] 
70 [dB] 


C. Case Definition 


A case code is comprised of one character and four digits, (for example C2461). 


It was used to define the basic MP assumptions. These include homogeneous versus 


IH propagation, the ocean conditions, the estimation resolution, The DR grid data 


and some additional parameters. The code is used in the computation of TDOA 


surfaces over the depth range (DR) grid. 
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The parameters are defined as follows. The leading character indicates the 
type of case. 

U - Homogeneous propagation, perfectly resolvable and identifiable MP. This 
is an idealized case where Eq. (1.2) is used to compute the unsorted (hence the U) 
time delays 7,72. It is used mainly for methodology purposes, demonstrating the 
differences between the TDOA and modified TDOA. 

S - Homogeneous propagation, perfectly resolvable but non identifiable paths. 
Eq. (1.2) is again used to generate the time delays which are then sorted according 
to delay magnitude (hence the S) at each DR point and denoted by ¢; and fg. 
The variable t; is now always the shortest, but it can no longer be associated with 
surface reflection nor can the longer t2 be associated with the bottom bounce. This 
case 1s useful as an idealized reference for the later IH limited resolution since they 
are conceptually similar. 

A,B,C,D - The realistic IH limited resolution non identifiable cases. The letter 


indicates the receiver time delay resolution in milliseconds as follows. 


A - 90.05 C - 0.5 
Bo =m D - 1.0 


The four digits define the parameters required for the numerical generation of 
the TDOA surfaces as follows. The first digit indicates the SVP. Four SVP’s of 


varying inhomogenuity are used throughout this research as shown below: 


case lst digit de a 3. 4. 
mae 
€ gs = C C (m/sec) 
lOO 
500 


gma, F=0.5 a? 0. a. 0.10$ 


The SVP digit is of course meaningless for the S and U cases where homuge- 
neous propagation (SVP digit = 1) is always assumed. 


The second digit details the observer and ocean depth information. 


TABLE E.1 
DEPTH DATA 


case 2nd digit 1 ji 3 4 ) 
Ocean depth m 503 Ore - 518 550 
Observer depth 162 162 - 162 162 


The third and fourth digits describe the discrete depth and range grid used 


for the direct function table. 


TABLE E.2 

GRIP DATA 
Grid density 1 2 3 4 9) 
Initial depth [m] 160 - 160 160 160 
Depth Step size [m] —50 25 —100 10 
Number of steps it 724 6 11 
Initial range [m] 500 500 500 500 500 
Range step size [m] 500 1000 500 1000 250 
number of steps 25 20 49 13 99 
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APPENDIX F 
THE EIGENRAY ACOUSTIC MODEL 


This appendix describes the selection of the eigenray model and the evaluation 
of the SMART model. 


A. REQUIREMENTS AND MODEL SELECTION 
The general goal of the eigenray model is to find and trace all the rays passing 
through the twe s.«cified receiver and source end points. The following is required 
from an eigenray model if it is to be used for MP tracking. 
1. Correctness and High Accuracy 
The total sound travel time from the source to the receiver may be of 
the order of a few tens of seconds. The milliseconds of travel time difference, 
while vital for the MP ranging, are a very small percent of the total travel time. 
This is especially true for the longer ranges where the travel time grows while the 
differences decrease. The resulting output accuracy requirement is of the order of 
001% which dictates even higher precision in the model’s internal computations. 
The issue of correctness relates to the possibilities of not finding an existing path 
or producing imaginary nonexisting paths; it too 1s related to the computer word 
length. 
2. Compactness 
The eventual use of a MP tracking algorithm in a real time systems appl- 
cation imposes additional constraints on run time, supporting hardware, and size 
of the model. Research oriented programs requiring extra long computer words or 


computation time will not fit the smaller data processors usually available for 
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target tracking in onboard systems. What is required is a program running on a 


microcomputer that can solve a typical eigenray problem in a matter of seconds. 


B. MODEL SELECTION 

Some eigenray models were developed in recent years, but the conflicting re- 
quirements of accuracy and compactness, narrows the selection to two programs. 
One model called SMART, for SMall Acoustic Ray Tracer was developed by Novick 
of Mission Rain os Corp. [Ref. 22]; the other is a research tool developed by Spies- 
berger as an application of an earlier model [Ref. 31]. 

The SMART model was selected because it was better integrated as a working 
model and provided more detailed output data. The Spiesberger model has is the 
advantage that it can account for a sloped ocean bottom. However this model lacks 
a reliable transmission loss output and is not well integrated for the intended use 
here. 

The eigenray model is used in a fashion that is separable from the rest of the 
MP inversion or target tracking algorithm. This provides for easy replacement or 
upgrading of the model. This also led to our decision to use the more reliable 
integrated model for this first IH compensation attempt, and maintain the option 
to upgrade and extend it in the future. 

The SMART model was evaluated extensively and the main result of this 


evaluation is presented in Section D. 


C. SMART MODEL - A BRIEF DESCRIPTION 
Finding the eigenrays is the first task accomplished by the model. Source and 


receiver depths are set to the specified input parameters and an iterative search 
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for the initial angle of the ray which will reach the receiver depth at the specified 
range, is conducted. The interval of possible ray angles is divided into sectors for 
which the dependence of range on the angle is a continuous well-behaved monotonic 
function. 

Figenrays are searched for in every sector that spans the desired range by 
successive approximation. the procedure terminates when the eigenray reaches a 
specified range accuracy. Only rays with a total number of bounces less than a 
specified input value are considered. Forward ray tracing equations similar to those 
presented in Chapter Three are used to locate the limits of the angle sector, to try 
out angles in the search process, and to compute the detailed output for every 


eigenray found. 


D. LIMITATION OF THE MODEL 

The program does not find eigenrays that are within +0.15 deg from the 
horizontal. This is partially caused by the limitation of intensity computation at 
these angles (recall discussion in Chapter 3 and Eq. (3.20)). Another cause for 
the limit is the finite computer word length and the increment used for searching 
the ray around the horizontal direction. The +0.15° limitation relates to both 
angle at the source and at the receiver. An example of the problem is shown in 
Fig F.1. The bounce count (BC) of case C1111* computed using SMART shows 
a high value along the receiver depth at long range. No loss of a path or loss of 


* 


resolution should give rise to this phenomena.** The unexpected BC was traced 


* See Appendix E for case definition. 
** Near the surface and the bottom the corresponding single bounce is not 
resolvable from the direct path. At depth of —170 m the single surface and 
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of the eigenray model to detect the direct path. This results from the direct ray 
being close to horizontal at those ranges. (The vertical separation between the 
grid depth line at 152 m and receiver at 162 m is 10 m.) These limitations were 
brought to the attention of the author of the SMART model and are currently 
being corrected. Since the model is employed in the prefilter in a modular way, 
these limitations do not effect the generality of the proposed inversion procedure. 

The choice of The SVP and geometry selected for use in this work was de- 
signed to avoid this problem by limiting the grid to a region were the rays are not 
horizontal at all. This was done by (1) restricting the DR grid to +60 - +160 m 
while the receiver was placed at 162 m; (2) limit the maximum range to 25 km; 
and (3) restricts the SVP gradients to positive surface gradients only. The overall 


geometry is shown in Fig. F.2. 


bottom bounce are not resolvable. 
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BOUNCE COUNT 
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Fig. F.1. Bounce count with missing horizontal rays. 
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Fig. F.2. Grid without horizontal rays. 
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APPENDIX G 
THE LOCAL SEARCH ALGORITHM 


In this Appendix the component of the local search algorithm, the inversion 
test, the inward search and the outward search, are described. A pseudo code 
description of'the algorithm is also included. The input and output parameters 
used by the algorithm are defined in Chapter Three. 

A Numerical analysis procedure, developed for finding zeros of complex func- 
tions [Ref. 29] was modified to apply to the search. For this purpose the search has 
to be viewed in the context of the graphical statement of the problem as discussed 
in Chapter Three, that is, as finding the DR coordinates of the point of intersection 


of the constant t, and t2 contours (see Fig. 3.12). The inversion problem 


(a)=() wy 
(x) (4) =(0) 


Where the D and R are the unknowns. If the vector function f is written as two 


can be rewritten as 


separate functions f; and f2 Eq. (G.2)can be written as the two equations 
D-f;' ae = 0 (G.3a) 
2 


The contour lines of Fig. 3.12 represent loci of the DR solution of Eq. (G.3). 


Every point p on the DR grid with depth D, and range R, can be defined in terms 


2a 


of it’s relation to the input t; contour by a one digit number C,; given by 
il stip >t 
Chi = ‘0 stip =h1 (G.4) 
and in terms of the it’s relation to the second input contour t2 by another digit 
C2 given by 
1 ; top > to 
Cy2 = ‘0 top = 22 (G.5) 
Cy, and Cy2 are referred to as characteristic digits. An example of the two char- 


acteristic digits for some points is shown on Fig. G.1. 


A. THE INVERSION TEST 

A 1x1 rectangle can contain a solution only if its characteristic digits Cp 
and C’,2 are either zero or if they acquire at least two different values in one of the 
four rectangle corners. 

Some examples are shown in Fig. G.1(c). In (a) both Cy; and Cy2 ave 1 and 
—1 at different points. In (b) C,; is 0 and 1, while C,2 is 1 and —1. This rule is a 
simplified version of a rule analytically derived in Ref [29] from Cauchy’s principal 
value argument. 

Note that this simple test is anecessary but not a sufficient condition. Fig. G.1c 
shows an example of a rectangle which succeeds the test but does not contain a 
solution. If this simple conditional test succeeds it 1s followed by a second more 
rigorous test which checks to see if the two contours actually intersect or just pass 
through the rectangle. This is done by actually attempting to locate the solution, 


by means of the inward search described later. 
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Fig. G.1. Examples of characteristic digits. 
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B. THE OUTWARD SEARCH 

The outward search is a recursive procedure. It performs the inward search 
over successively growing search rectangles. If a call to the inward search returns 
in a failure (no solutions in the current rectangle) the search area is increased by 
tripling the sides of the current rectangle. A new current rectangle is thus formed 
centered around the previous one with an area nine times larger. The outward 
search is now called recursively to search the new rectangle. Again it tests for 
existence of a solution by calling the inward search. Note that in its growth the 
outward search may in one step include two or more solutions which would then 
all be located by the inward search. 

The search rectangle area is increased by a 9:1 ratio by means of tripling the 
length of both its sides. This is done in order to minimize the duplicated effort 
associated with re- examining the original rectangle as part of the recursive search 
of the new rectangle. However the rectangle is increased only up to the portion 
of the maximum search area size that still remains within the limits of the grid. 
If the inward search for the current rectangle is unsuccessful and the search area 
cannot be further increased in any direction, the procedure terminates in a failure, 
without any solution. 

Upon termination the search procedure outputs the size of the current area, 
the number of solutions found and the solution quad triangles. If more than one 
solution is found the calling inversion procedure ~vi!l select the most likely one, i.e., 
the one closer to the predicted position after the interpolation is performed. An 
example of the growing search rectangles is shown on Fig. 3.18 in the dashed line. 
The search is initiated at point D,, Rp, and the existence of solution is detected 


after two iterations. 
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C. THE INWARD SEARCH 

The inward search attempts to locate a solution in a given input rectangular 
search area. It starts by performing the conditional inversion test on the complete 
area. If the test succeeds, the region is divided into four (quadrant) rectangles 
by halving both of its axes. Each quadrant is then searched for solutions by a 
recursive call to the inward search*. If the possible existence of a solution in the 
quadrant is indicated by the inversion test the quadrant is subdivided and the 
process continues recursively. If the test of the quadrant fails, the search in it 
terminates. The inward recursion thus terminates when either the search area 
reaches the minimum rectangle (of size 1 x 1), or if no solution exists in the a 
specific branch (of the B tree with branching factor 4). 

The triangle(s) in the 1 x 1 rectangle is (are) found by dividing the rectangle 
into it’s two DR grid triangles and performing a test on each to see if it surrounds 
the point. The surrounding test uses vector product to verify that the input t,, te 
point is on the same side of all lines that make the sides of the triangle. The test is 
as follows: Three vectors S;, S_ and S3 are defined in the t;t2 plane as the three 
sides of the triangle in a cyclic order (see Fig. G.2). Three corresponding vectors 
V,, V2 and V3 are formed by the input t; and ¢2 point and the three angles of the 
triangle A,,A2,A3 respectively. The side of the input point relative to a “side 
vector” S is given by the sign of the vector product of the corresponding vectors 
S and V vectors. An input point inside a triangle will be on the same side of all 


the three side vectors. That is 


sign [Vi x Si] = sign [V2 x S| = sign [V3 : S3] (G.6) 


* The search is therefore a B tree search with a branching factor of 4. 
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Fig. G.2. Vectors for inside test. 
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In Fig. G.2 the six vectors are plotted for a point inside (p, and full line) and for 
a point outside (pg and dashed lines) a triangle. The signs marked along the V 
vectors are the signs of the corresponding vector product (V; x S;). 

The solutions from all the branches of the recursion are reported by the inward 
search. This ensures that all solutions found in the searched area centered around 
the reference point are reported for further consideration. In the example shown 
in Fig. 3.18 three iterations of the inward search are required to locate the actual 


triangles that surround the point. In this case there were two solutions. 


Pseudo code of the local search algorithm follows. 
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D. PSEUDO CODE FOR THE LSA 


TYPE box : a square area in the depth range grid defined by 


the DR indexes. 


FUNCTION inside ( tl, t2 ,tri): boolean; 
inside :- 1 if the point t1,t2 is inside the triangle defined 
by the three quads of tri on the t1t2 plane; 


O otherewize. 


FUNCTION close ( INVAR: depth, range): box 
close :- the 1 j (depth range ) indexes of the quad bex::on 
the depth range grid surrounding the input 


depth range point. 


ENDFUNCTION; 
PROCEDURE inversion test ( INVAR : t1, t2, current; 
OUTVAR : result ) 
IF 
(CP1l equals O or changes sign in the current area) 
AND (CP2 equals O or changes sign in the current area) 


THEN result :-1 ELSE result :- 0; 


PROCEDURE owtward ( INVAR : t1, t2, current , maxsize; 
OUTVAR :Qlist, number_of_solutions, areasize) 
BEGIN 
IF current = previous THEN exit ;reached the size limit. 
ELSE 
BEGIN 


inward ( tl, t2,current,glist, num of solutions) ; 


Zee 


IF number _of_ solutions =/= 0 THEN exit;Solution found. 
Else 
BEGIN ; Recourse. 
current :- triplsized (current). 


outward(t1,t2,current,num_of_ solution, areasize) 


END; 
ENDIF; 
END; 
ENDIF; 
ENDPROCEDURE. 
PROCEDURE inward (INVAR ¢; t1,t2,current 
OUTVAR : qlist,num_of_solutions) ; 

BEGIN 


IF current = 1 x 1 THEN 
split current box into triangles tril, tri2; 
FOR tri IN (tril, tri2) DO 
IF inside (t1,t2,tri) THEN 
add tri to qlist; 
num_of_solutions:- num_of_solutions+1; 
ENDIF; 
ENDFOR; 
ELSE 
qlist :- empty; 
num_of solution :- 0; 
IF not(inversion test(tl, t2, current)) THEN exit 


ELSE 


split current into four quadrants 1 thru 4; 
FOR next IN quadrants 1 through 4 DO 
inward ( t1,t2,next,list,num); 


qlist :- glist + list; 


num_of_ solution :- num_of_solution + 1; 
ENDFOR; 
ENDIF; 
ENDIF; 


ENDPROCEDURE. 


PROCEDURE inversion ( INVAR :t1, t2, DRref,maxsize 
OUTVAR : depth, range,aver_num,aver_areasize) 
var current : box; 
solutionlist : list of depth range solutions; 
BEGIN 
current :- close( DR refpoint) ; 
outward ( t1,t2,current,qlist,num_of_solutions,areasize) ; 
IF num_of_solution = O THEN depth, range :- DRref ;a failure 
ELSE 
BEGIN 
aver num :- average_nun x a + num_of solutions x (1-a); 
aver_areasize :- aver _areasize x a + areasize x (1 - a); 
FOR tri IN qlist DO 
interpolate depth range using tri; 
append solution to solution list; 
ENDFOR; 
ENDIF; 
select depth range solution closest to DRref point as depth 


range values to output. 


ENDPROCEDURE. 
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